A0-/U12  237  OHIO  STATE  UNIV  COCUNRUS  DEPT  Of  GEODETIC  SCIENCE  A--ETC  F/G  A/* 

THE  DO»N»AAO  CONTINUATION  TO  THE  EARTH'S  SURFACE  OF  TRUNCATED  S-- ETCIHI 
DEC  61  C  JEKELI  Fl962A-7q-c-0rt?7 

UMCLASSIF1EO  323  AFSL-TR-61-0361  ML 


MICROCOPY  RESOLUTION  TEST  CHART 

NATIONAL  BURIAU  OF  STANDARDS  1<>6J  A 


ADA  11223  7 


\ 

1 

1  AFGL-TR- 81-0361 


THE  DOWNWARD  CONTINUATION  TO  THE  EARTH'S  SURFACE  OF  TRUNCATED 
SPHERICAL  AND  ELLIPSOIDAL  HARMONIC  SERIES  OF  THE  GRAVITY  AND 
HEIGHT  ANOMALIES 


Christopher  Jekeli 


The  Ohio  State  University 
Research  Foundation 
1958  Neil  Avenue 
Columbus,  Ohio  43210 


December,  1981 
Scientific  Report  No.  11 

Approved  for  public  release;  distribution  unlimited 


r  AIR  FORCE  GEOPHYSICS  LABORATORY 

AIR  FORCE  SYSTiUS  COMMAND 
*  UNITED  STATES  AIR  FORCE 

HANSCOM  AFB,  MASSACHUSETTS  01731 


gg  08  .  ol7 


UNCLASSIFIED 


SECURITY  CLASSIFICATION  of  THIS  PAGE  (When  Dele  Entered ) 


REPORT  DOCUMENTATION  PAGE 

READ  INSTRUCTIONS 

BEFORE  COMPLETING  FORM 

1.  REPORT  NUMBER  2.  GOVT  ACCESSION  NO. 

AFGL-TR-8 1-0361  -  A  VA 

3.  RECIPIENT'S  CATALOG  NUMBER 

«.  TITLE  (end  Subtitle) 

THE  DOWNWARD  CONTINUATION  TO  THE  EARTH'S 
SURFACE  OF  TRUNCATED  SPHERICAL  AND  ELLIP¬ 
SOIDAL  HARMONIC  SERIES  OF  THE  GRAVITY  AND 
HEIGHT  ANOMALIES 

S.  TYPE  OF  REPORT  A  PERIOD  COVERED 

Scientific  Report  No.li 

6.  performing  O^G.  REPORT  number 

Dept  .of  Geodetic  Sci.  #323 

7.  AUTHORfaj 

Christopher  Jekeli 

8.  CONTRACT  or  grant  NUMBERfs) 

F19628-79-C-0027 

9.  PERFORMING  ORGANIZATION  NAME  AND  ADDRESS 

The  Ohio  State  University 

Research  Foundation 

1958  Neil  Ave.,  Columbus,  Ohio  43210 

10.  PROGRAM  ELEMENT.  PROJECT,  TASK 

AREA  &  WORK  UNIT  NUMBERS 

661102F 

2309G1AW 

It.  CONTROLLING  OFFICE  NAME  AND  ADDRESS 

Air  Force  Geophysics  Laboratory 

Hanscom  AFB,  Massachusetts  01730 

Contract  Monitor:  George  Hadgigeorge  LWG 

12.  REPORT  DATE 

December,  1981 

13.  NUMBER  OF  PAGES 

140 

14.  MONITORING  AGENCY  NAME  A  ADDRESS (II  dllterent  tram  Controlling  Ottice) 

IS.  SECURITY  CLASS,  (ol  thie  report) 

Unclassified 

15a.  DECLASSIFICATION/ DOWNGRADING 
SCHEDULE 

16.  DISTRIBUTION  STATEMENT  (cl  thte  Report) 

A  -  Approved  for  public  release;  distribution  unlimited 

17.  DISTRIBUTION  STATEMENT  (ol  the  ebetrect  entered  In  Bloch  30,  If  different  from  Report) 

IS.  SUPPLEMENTARY  notes 

IS.  KEY  WORDS  (Continue  on  rovoreo  tide  if  necoeeery  and  Identify  by  block  number) 

geodesy,  spherical  harmonics,  gravity,  potential  theory 

20.  ABSTRACT  (Continue  on  reveree  tide  It  neeeeeary  and  Identify  by  block  number) 

The  problem  of  divergence  of  the  geopotential  spherical  harmon¬ 
ic  series  at  the  earth's  surface  is  investigated  from  a  numerical, 
rather  than  a  theoretical  approach.  It  is  shown  that  previous  num¬ 
erical  evaluations,  based  on  the  expansion  of  the  potential  gener¬ 
ated  by  the  fflasses  between  the  point  of  computation  and  the  bounding 
sphere,  are  inconclusive  with  respect  to  the  magnitude  of  the  down¬ 
ward  continuation  error  that  is  attributable  to  series  divergence. 

A  more  representative  model  of  the  earth's  potential  is  devised  on 

DD  I  jan *7i  1473  edition  of  1  NOV  ss  is  obsolete  UNCLASSIFIED 


SECURITY  CLASSIFICATION  OF  THIS  PACE  flWiMi  Dete  Entered) 


I 


UlYCi-.jiOOi.r  lau 

SECURITY  CLASSIFICATION  OF  THIS  RAOEflWin  Dsn  Bnfrmd)  _ 

the  basis  of  a  density  layer,  which,  in  the  spherical  approximation, 
generates  a  gravity  field  whose  harmonic  constituents  decay  accord¬ 
ing  to  an  accepted  degree  variance  model.  This  field,  expanded 
to  degree  300,  and  a  topographic  surface  specified  to  a  corres¬ 
ponding  resolution  of  67  km  are  used  to  compute  the  differences 
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anomalies  at  the  surface  of  the  earth  model.  Up  to  degree  300, 
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1 .  Introduction 


Recent  years  have  seen  a  guantua  juap  in  the  continuing 
efforts  to  improve  and  expand  our  knowledge  of  the  earth's 
gravity  field.  It  is  Barked  by  the  transition  froa 
terrestrial  aeasureaents,  United  essentially  to  land  areas, 
to  aeasureaents  in  outer  space  which  are  used  to  ascertain 
the  global  gravity  field.  The  foremost  instrument  advancing 
this  effort  has  been  the  satellite  borne  altiaeter,  enabling 
a  direct  aeasureaent  of  an  equipotential  surface,  the  geoid, 
over  such  of  the  oceanic  surface  of  the  earth.  Further 
strides  will  undoubtedly  be  made  by  the  planned  GHAVSAT 
mission  (see  the  report  by  the  National  Research  Council, 
1979)  utilizing  the  measurements  of  sa tellite- to-satellite 
tracking  and  achieving  coverage  over  the  entire  globe. 
Looking  ahead  into  the  not  so  distant  future,  satellite 
borne  gradiometers  will  provide  even  greater  detail  and 
accuracy.  While  the  obvious  advantage  of  measuring  the 
earth's  gravity  field  at  satellite  altitudes  is  global 
accessibility  within  a  relatively  short  period  of  tiae,  the 
fundamental  difficulty  is  the  translation  or  "downward 
continuation"  of  the  data  to  the  earth's  surface  where  they 
are  most  needed.  In  principle,  several  procedures  to 
achieve  this  translation  exist;  all  rely  to  some  extent  on  a 
simplifying  assumption  such  as  a  perfectly  spherical,  or  a 
flat,  earth.  On  account  of  the  enormous  amount  of  data  that 
satellite  missions  provide,  the  number  of  methods  to 
simultaneously  process  the  entire  data  set  is  reduced 
considerably.  The  method  of  harmonic  analysis  of  the 
gravitational  potential  will  come  under  close  examination  in 
this  paper.  Its  feasibility  froa  the  coaputational 
standpoint  cannot  be  easily  challenged,  even  for  extremely 
dense  data  coverage.  However,  far  froa  being  a  panacea,  it 
is  also  associated  with  several  problems.  Aside  froa  an 
instability  in  the  propagation  of  noise,  the  aost  nagging 
question  is  one  of  correct  theory.  It  is  the  latter  which 
will  be  studied  here,  not  by  delving  into  areas  of  pure 
theory,  bot  rather  on  a  numerical  basis,  which,  it  is  felt, 
will  provide  soae  value  to  the  scientist  who  aust  eventually 
make  use  of  the  data. 
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1.1  The  Problem  and  Background 

Proa  classical  potential  theory,  we  know  that  the 
(Newtonian)  potential  due  to  attracting  Basses  is  an 
haraonic  function  in  free  space.  That  is,  its  second 
derivatives  are  continuous  and  it  satisfies  Laplace's 
equation;  moreover,  it  is  regular  at  infinity  (Gunter,  1967, 
p.25).  Kellogg  (1903,  p.220)  shows  that  an  haraonic 
function  is  also  analytic  in  its  region  of  haraonicity  (cf . 
Cauchy  Biemann  equations  in  the  theory  of  complex 
variables) .  The  solution  to  the  exterior  Dirichlet  boundary 
value  problem  states  that  given  values  of  the  potential 
everywhere  on  a  known  surface  enclosing  the  masses  it  is 
determined  uniguely  in  the  space  exterior  to  that  surface. 
Applied  to  earth  orbiting  satellites,  the  Known  surface  is 
the  sphere  that  contains  the  satellite  orbits.  Because  the 
potential  is  also  analytic  in  the  region  between  the  earth's 
surface  and  the  orbital  sphere  (we  remove  the  atmosphere, 
see  below) ,  by  the  uniqueness  of  analytic  continuation,  the 
potential  function  outside  the  orbital  sphere  represents  the 
potential  in  the  entire  region  above  the  earth's  surface, 
i.e.  in  the  largest  region  wherein  the  actual  potential  is 
analytic. 

Por  an  irregular  density  distribution  such  as  the 
earth's,  a  closed  form  of  the  potential  in  space  cannot  be 
found.  Instead,  it  is  often  represented  as  a  series,  in  one 
fora  or  another;  however,  any  series  is  associated  with  a 
particular  region  of  convergence  and  cannot  converge  to  the 
true  potential  in  the  total  space.  Since  we  are  dealing 
with  exterior  potentials,  the  region  of  convergence  is  an 
exterior  region  that  contains  the  point  of  infinity,  and  it 
is  separated  from  the  region  of  divergence  (the  interior 
region)  by  the  so-called  surface  of  convergence.  In 
general,  the  region  of  convergence  nay,  or  may  not,  contain 
the  maximum  region  of  analyticity  of  the  potential,  nor  is 
the  surface  of  convergence  necessarily  a  sphere.  These 
facts  were  convincingly  demonstrated  by  Krarup  (1969)  and 
Moritz  (197b);  see  below,  in  many  cases,  it  is  possible  to 
derive  a  series  which  converges  to  the  potential  in  the 
region  where  the  outer  series  diverges;  we  call  this  the 
inner  series. 

Owing  to  the  near  spherical  shape  of  the  terrestrial 
body,  the  most  familiar  series  is  the  spherical  harmonic 
series.  Kellogg  (1953,  p.193)  showed  that  the  spherical 
harmonic  series  converges  uniformly  to  the  potential  outside 
any  sphere  containing  all  the  attracting  masses  and  centered 
at  the  origin  of  the  coordinate  system.  Strictly,  this 
theorem  finds  no  application  in  our  physical  world  since  it 
guarantees  convergence  only  outside  the  sphere  enclosing  the 
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entice  universe.  Within  our  limited  scope  o f  terrestrial 
and  near-earth  applications,  however,  the  masses  outside  our 
solar  systea  have  negligible  gravitational  effect. 

Moreover,  we  may  si* ply  redefine  the  exterior  gravitational 
field  of  the  earth  with  appropriate  corrections  so  as  to 
exclude  the  effects  of  an  atmosphere  and,  say,  of  the  sun 
and  the  moon .  The  latter  and  possibly  other 
extra-terrestrial  bodies  are  mathematically  aoved  to 
infinity  (where  they  have  no  gravitational  influence)  by 
subtracting  the  corresponding  tidal  potential.  The 
ataosphere  is  oust  conveniently  embedded  (conceptually) 
inside  the  earth  such  that  the  center  of  *ass  remains 
undisturbed  (Moritz,  1974) ;  of  course,  the  resulting  change 
in  the  earth’s  exterior  gravity  field  *ust  be  accounted  for 
when  coaparing  terrestrial  data  with  data  downward  continued 
froa  satellite  altitudes.  Positioning  the  origin  of  our 
coordinate  systea  approximately  at  the'  earth’s  center  of 
■ass,  we  therefore  have  guaranteed  convergence  of  the 
spherical  harmonic  series  of  the  potential  outside  the 
sphere  whose  radius  equals  the  farthest  distance  of  the 
earth’s  surface  from  the  earth’s  center;  this  is  the  top  of 
the  Chimborazo  fountain,  in  central  Ecuador  (latitude  -1?4) , 
with  a  radial  distance  of  about  6384403  a  (Sjoberg,  1977). 

A  aore  general  result  was  rigorously  proved  by  Krarup 
(1969,  chapter  3),  naaely  that  the  spherical  haraonic 
expansion  converges  everywhere  on  and  outside  the  smallest 
sphere  (called  the  limit  sphere)  that  contains  all 
singularities  of  the  potential  and  its  analytic 
continuation.  Thus,  the  proof  that  the  potential  series 
converges  everywhere  at  the  earth's  surface  would  be 
complete  if  the  potential  could  be  analytically  continued 
down  to  the  so-called  Bjerhaaaar  sphere  (the  sphere  that  is 
entirely  enclosed  within  the  earth).  Kellogg  (1953,  p.197) 
comments  that  the  potential  function  representing  the 
potential  of  an  analytic  density  distribution  bounded  by  an 
analytic  surface  can  be  oontinued  analytically  across  the 
surface.  Of  course,  in  view  of  Poisson’s  equation,  the 
actual  potential  is  not  represented  by  this  function  at 
points  of  nonzero  density;  indeed,  its  discontinuous  second 
derivatives  preclude  its  being  analytic  on  the  surface.  The 
possibility  of  analytically  continuing  the  geopotential 
inside  the  irregular  masses  of  the  earth  seems  very 
doubtful,  for  as  Krarup  points  out,  if  it  were  possible  for 
some  given  mass  distribution,  the  mere  addition  of  a  mass 
point  ("grain  of  sand")  above  the  limit  sphere  introduces  a 
singularity  in  the  potential  function  at  this  point  and 
thereby  destroys  the  analyticity  of  the  continuation. 
Therefore,  given  that  the  series  for  the  potential  converges 
with  certainty  only  outside  the  bounding  sphere,  the 
question  arises  whether  there  exists  any  justification  for 
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usiou  the  series  a t  or  near  the  surface  or  the  earth. 

Claias  of  both  proof  and  disproof  of  series  convergence 
at  the  earth's  snrface  appear  in  the  geodetic  literature; 
none  is  logically  sound.  The  proof  by  Arnold  (1978)  of 
convergence  everywhere  on  the  snrface  is  patently  flawed,  as 
is  Horrison's  (1970)  conjecture  of  divergence  everywhere 
(see  Appendix  A).  The  "grain  of  sand"  example  was  used  by 
noritz  (1980,  p.64)  to  argue  on  the  instability  of  the 
property  of  convergence,  iaplying  that,  in  the  strictest 
sense,  the  size  and  shape  of  the  surface  of  convergence,  if 
it  is  not  the  bounding  sphere,  is  not  well  defined.  Any 
further  theoretical  advances  on  the  behavior  of  spherical 
haraonic  series  near  the  surface  of  an  attracting  body  will 
cone  by  studying  density  distributions  bounded  by  surfaces, 
both  of  which  are  mathematically  regular  by  some  aeasure.  A 
study  which  approaches  this  type  of  analysis  is  the  one  by 
Kbolshevnitov  (1977) ,  who  finds  upper  bounds  for  the  decay 
rate  of  spherical  aaraonics,  on  the  bounding  sphere,  with 
respect  to  the  degree  n  for  variously  structured  bodies. 
These  upper  bounds  are  generally  proportional  to  (fixed) 
neqative  powers  of  n ,  depending  on  the  aeasure  of  regularity 
of  the  density  and  bounding  surface.  Such  decay  rates  are 
insufficiently  strong  to  yield  convergence  below  the 
bounding  sphere;  yet  as  they  are  only  upper  bounds, 
convergence  cannot  hereby  be  excluded. 

While  the  precise  convergence  surface  for  series 
corresponding  to  arbitrary  density  distributions  bounded  by 
arbitrary  surfaces  has  eluded  theorists,  several 
fundamental,  as  well  as  interesting,  results  have  been 
established.  Star up  (1969)  examined  the  potential  resulting 
froa  a  Kelvin  transformation  of  the  potential  of  a  uniform 
aass  distribution  on  a  straight  line.  Using  this  example, 
he  disproved  the  intuitive  notion  that  the  surface  which 
separates  regions  of  convergence  and  divergence  is  always  a 
sphere,  a  general  theory  regarding  the  shape  of  the  surface 
of  convergence  for  special  spherical  haraonic  series  was 
developed  by  Ecker  (1972).  Ke  proved  that  a  sphere  is  the 
surface  of  convergence  tor  rotationally  symmetric  potentials 
(i.e.  series  of  zonal  harmonics  only),  while  a  torus  defines 
the  surface  of  convergence  for  a  series  of  only  tesseral 
harmonics  (Krarup's  example).  Other  surfaces  of  convergence 
lying  between  these  two  extremes  result  for  series  of  only 
those  sectorials  whose  degree  n  and  order  a  satisfy  the 
relation  a=ln,  i  being  predefined  and  0i(  il.  The  case 
l  *0  represents  the  series  of  zonals,  and  Ecker  proved  the 
following  result: 
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Theoieui:  A  serits  of  spneiical  Harmonics  £  (  — )n  a  P  (cos6) 

n=0  r  n  n 

c ocvei'OtiS  everywhere  outside  the  sphere  of  radius 

a  .  where  ,  . 

r  lim  n  - r 

P  =  sup  /|ani 

and  diverges  almost  everywhere  inside  this  sphere. 


Here,  r,S  are  the  polar  coordinates,  radius  and  colatitude, 
and  denotes  the  Legendre  polynomial  of  n-th  degree.  The 
qualification  "almost  everywhere"  for  divergence  is  included 
to  allow  convergence  on  sets  of  aeasurt  zero  (regions  having 
no  volume)  inside  the  convergence  sphere  (e.g.  any  series  of 
odd  zonals  converges  for  every  r > 0  at  the  equator,  since 
Pi—,<0>=0)  -  The  radius  of  convergence  for  the  by  now 
classic  example  of  a  homogeneous  oblate  ellipsoid  of 
revolution  (Jung,  1956,  p.54  3;  floritz,  19B0,  p.52)  is  found, 
using  the  above  theorem,  to  be  f  =B,  the  focal  distance  of 
the  ellipsoid,  showing  also  that  the  surface  of  convergence 
may  not  bound  the  generating  masses  (i.e.*  in  this  case  the 
potential  function  can  be  anaytically  continued  to  the 
sphere  of  radius  B  but,  of  course,  does  not  represent  the 
potential  inside  the  ellipsoid) . 

The  question  of  convergence  or  divergence  of  the 
potential  at  the  earth's  surface  may  be  circumvented  by  the 
Hunge-Krarup  theorem  (Krarup,  1969;  Horitz,  1980,  p.67) . 
Briefly,  this  theorem,  already  known  to  Walsh  (1929,  p-535) 
for  the  inner  potential,  states  that  a  function  harmonic 
outside  the  earth's  surface  may  be  approximated  arbitrarily 
well  in  its  region  of  harmonicity  by  a  function  which  is 
harmonic  outside  a  given  sphere  totally  inside  the  earth. 
Obviously,  the  spherical  harmonic  series  of  an  harmonic 
function  converges  everywhere  outside  any  sphere  contained 
entirely  in  its  region  of  harmonicity,  in  particular  on  the 
earth's  surface  if  this  sphere  is  embedded  entirely  within 
the  earth.  It  should  be  noted  that  the  Runge-Krarup  theorem 
is  an  existence  theorem;  it  guarantees  only  the  existence  of 
an  approximating  function  and  does  not  provide  the  method  to 
find  it-  Furthermore,  nothing  is  said  about  the  closeness 
of  corresponding  individual  terms  of  the  two  series  for  the 
actual  and  approximating  potentials.  The  approximation  is 
arbitrarily  accurate  only  in  the  limit,  i.e.  for  the  total 
sum  (however,  one  can  expect  that,  because  of  the  near 
sphericity  of  the  earth's  surface,  the  corresponding  lower 
degree  terms  of  the  two  series  do  not  differ  substantially) . 
We  may  also  note  that  the  theorem  holds  for  any  exterior 
potential  no  matter  how  badly  its  series  diverges  below  the 
bounding  sphere,  mo  that  the  application  of  the  theorem  is 
not  contingent  on  the  instability  of  the  convergence  surface 
of  the  actual  potential  series. 
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The  fact  that  the  Runge-Krarup  theorea  says  nothing 
about  the  accuracy  of  individual  terns  of  the  potential 
series  is  very  restrictive  on  the  possible  use  of  the 
theorea  to  the  practicing  geodesist,  hence,  for  example, 
Sjoberg's  (1979)  statenent  that  "the  coefficients  of  the 
approximating  potential  ...  can  be  selected  arbitrarily 
close  to  the  coefficients  of  the  external  potential  ...  to 
any  desired  degree"  is  somewhat  misleading  if  not 
inaccurate.  For  tne  closer  the  coefficients  of  the 
approximating  potential  are  chosen  to  those  of  the  actual 
potential,  the  further  is  the  postponement  of  the  arbitrary 
closeness  of  the  approximating  series  to  the  true  potential. 
However,  for  the  example  of  a  point  mass  situated  on  the 
equator  of  an  oblate  spheroid  (so  that  the  ensuing  potential 
series  diverges  in  the  polar  regions),  Sjoberg  shows  that  in 
this  case,  the  postponement  is  not  unduly  exacerbated  if  the 
first  300  or  400  terms  o£  the  actual  divergent  series  are 
used  for  the  corresponding  terns  of  the  approximating 
potential,  aontz's  statement  (horitz,  1980,  p.66)  based  on 
the  Runge-Krarup  theorem  that  the  earth's  potential,  for 
practical  purposes,  can  always  be  considered  as  a 
"convergent  potential"  is  similarly  misleading.  This 
statement  should  not  be  interpreted  as  claiming  convergence, 
for  all  practical  purposes,  of  the  actual  series  of  the 
potential  on  and  outside  the  earth's  surface.  Instead  it  is 
a  statement  on  the  practical  equivalence  of  the  potential 
and  an  approximating  series  that  converges  on  and  above  the 
surface.  Term  lor  term,  especially  at  high  degrees,  the 
approximating  series  and  the  actual  series  must  be  quite 
different  since  the  former  converges  while  the  latter 
possibly  diverges  near  the  earth's  surface.  Therefore,  the 
Runge-Krarup  theorem  can  not  be  invoked  to  justify  the  use 
of  a  portion  of  the  series  of  the  actual  potential  at  the 
earth's  surface. 

Until  now  the  discussion  has  centered  on  the  convergence 
of  a  spherical  harmonic  series,  where  the  use  of  spherical 
coordinates  is  motivated  by  the  near  spherical  shape  of  the 
earth's  surface.  But  to  a  second  approximation,  the  surface 
of  the  earth  is  an  ellipsoid,  or  more  precisely  an  oblate 
spheroid,  an  ellipsoid  whose  equatorial  axes  are  equal  (i.e. 
it  is  a  surface  of  revolution)  and  whose  poles  are 
flattened.  The  question  arises  whether  the  use  of  different 
coordinates  such  as  ellipsoidal  coordinates  has  a 
significant  bearing  on  the  problem  of  convergence  at  the 
earth's  surface.  The  ellipsoidal  coordinates  for  which  the 
general  triaxial  ellipsoid  is  a  coordinate  surface  (a 
surface  defined  by  the  fixed  value  of  one  coordinate,  in 
tnis  case  one  or  the  semi-axes)  are  rather  more  difficult  to 
work  with  than  spherical  coordinates,  but  expansions  of  the 
potential  and  the  gravity  anomaly  have  been  formulated  in 
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terms  of  the  orthogonal  Lane  functions,  see  (Hobson,  1965, 
pp  473-475;  Walter,  197C;  and  Savrov,  1974).  Because  the 
deviations  between  best  fitting  triaxial  and  biaxial 
ellipsoids  are  of  the  sane  order  of  aagnitude  as  geoid 
undulations,  the  triaxial  ellipsoid  has  generally  been 
abandoned  as  an  approximation  to  the  earth *s  surface.  Thus, 
the  coordinate  system  to  be  considered  is  the  spheroidal 
system  (Hobson,  1965,  p.421,  see  also  chapter  3)  in  which 
the  expansion  of  the  potential  is  in  terns  of  faailiar 
Legendre  functions.  Though  not  widely  used  in  practice, 
these  coordinates  have  received  considerable  attention  in 
geodesy,  in  particular,  by  Jung  (1956)  and  Hotine  (1969)  . 

The  tern  "spheroid"  in  the  geodetic  volcabulary 
conventionally  denotes  an  eguipotential  surface  of  sone 
nornal  (reference)  potential.  On  the  other  hand,  the  tern 
"ellipsoid"  usually  implies  oblate  spheroid;  hence  it  will 
also  be  used  here  to  mean  exclusively  an  ellipsoid  of 
revolution  flattened  at  the  poles. 

Iaposing  rotational  symmetry  with  respect  to  one  of  the 
ellipsoidal  coordinates,  namely  the  longitude,  in  this  case 
yields  ellipsoidal  harmonic  functions  whose  structure 
differs  from  their  spherical  counterparts  only  in  the 
dependence  on  the  distance  from  the  origin.  Due  to  the 
corresponding  similarity  to  spherical  harmonic  expansions  we 
have  theorems,  such  as, 

1.  the  ellipsoidal  harmonic  series  of  the  potential 
converges  uniformly  everywhere  outside  the  ellipsoid  that 
bounds  the  generating  masses;  and 

2.  there  exists  an  ellipsoidal  harmonic  series  which 
converges  uniformly  above  the  "Bjerbaaoar  ellipsoid"  and 
approximates  the  potential  outside  the  earth's  surface 
with  arbitrary  accuracy. 

The  proof  of  the  first  statement  follows  immediately  from 
the  uniform  convergence  of  the  ellipsoidal  series  for  the 
reciprocal  distance  (equation  (3.27);  cf.  Kellogg,  1953, 
p.143,  and  see  also  Hobson,  1965,  pp. 430-433).  The  second 
statement  is  merely  a  corollary  to  the  Bunge-Krarup  theorem 
since  the  Bjerhammar  sphere  enters  only  to  relate  the 
theorem  to  spherical  harmonic  expansions  and  can  easily  be 
replaced  by  the  Bjerhamaar  ellipsoid.  Other  corresponding 
theorems  with  respect  to  the  "limit  ellipsoid"  or  to  the 
convergence  surface  of  a  series  of  ellipsoidal  zonals 
undoubtedly  exist,  but  may  require  more  exacting  proofs. 

Unfortunately,  the  transition  to  ellipsoidal  coordinates 
does  not  solve  the  problem  of  series  convergence  at  the 
earth's  surface  since  its  distance  from  the  bounding 
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ellipsoid  can  still  be  6  to  7  ka  (for  the  bounding  sphere  it 
is  as  such  as  25  ka) .  However,  we  recognize  that  the 
convergence  problea  is  a  manifestation  of  the  choice  of 
coordinate  sjstea,  for  the  ellipsoidal  series  will  converge 
to  the  potential  in  regions  where  convergence  of  the 
spherical  series  is  doobtfol  (e.g.  the  polar  areas).  Also, 
oe cause  the  equatorial  radius  of  the  bounding  ellipsoid 
could  exceed  the  radius  of  the  bounding  sphere,  the 
spherical  series  will  converge  in  regions  where  the 
ellipsoidal  series  way  not,  see  Pig.  1.  The  dependence  of 


Figure  1:  Bounding  sphere  versus  bounding  ellipsoid. 
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the  region  of  convergence  on  tie  coordinate  systea  is  even 
aore  directly  illustrated  by  siaply  changing  the  coordinate 
origin,  at  which  the  bounding  sphere  is  (always)  centered. 
Therefore,  the  saallest  bounding  sphere  is  obtained  if  the 
origin  of  the  systea  of  spherical  coordinates  coincides  with 
the  geoaetrical  center  of  the  earth.  The  feasibility  of  a 
choice  of  coordinates  other  than  ellipsoidal  coordinates 
that  guarantees  convergence  of  the  corresponding  series 
significantly  closer  to  the  earth's  surface  seeas  unlikely 
since  the  next  approxiaation  to  the  surface  is  a 
considerably  aore  coaplex  geoaetric  figure.  This  is  an 
upward  continuation  of  the  telluroid,  or  soae  smoothed 
version  oi  it.  The  telluroid  (Heiskanen  and  Moritz,  1967, 
p.292)  is  the  surface  of  points  at  which  the  noraal 
potential  equals  the  gravity  potential  at  the  corresponding 
points  on  the  earth's  surface,  where  correspondence  is 
established  if  the  telluroid  and  surface  points  lie  on  the 
sane  (noraal)  pluab  line  in  the  noraal  gravity  field.  The 
telluroid  iaitates  the  earth's  surface  guite  closely  since 
their  difference,  the  height  anoaaly,  varies  as  smooth ly  as 
the  geoid  undulation  with  average  values  of  30  a.  However, 
the  corresponding  coordinate  systea  (in  which  the  bounding 
telluroid  is  a  coordinate  surface)  will  be  too  abstruse  to 
work  with. 

The  essence  of  this  paper  addresses  the  question  of 
whether  the  probable  divergent  character,  at  the  earth's 
surface,  of  the  spherical  harnonic  expansion  of  the 
disturbing  potential  (and  gravity  anoaaly)  eliainates  it 
froa  the  repertoire  of  viable  aethods  of  downward 
continuation.  In  light  of  the  foregoing  suaaary  of  the 
theoretical  viewpoints  on  convergence  and  divergence  of  the 
earth's  potential  series,  the  analysis  will  be  based  on  the 
assuaption  that  the  series  definitely  diverges  below  the 
bounding  sphere.  Although  divergence  has  not  been  proved, 
this  premise  is  certainly  reasonable,  if  only  as  the 
worst-case  situation . 


1.2  Other  Problems  and  Methods  of  Downward  Continuation 

The  question  of  divergence  or  convergence  at  the  earth's 
surface  will  never  impede  our  computational  abilities  in 
practical  situations.  That  is,  with  a  finite  nuaber  of 
aeasureaents  of  the  potential  we  can  deteraine  only  a  finite 
nuaber  of  coefficients  of  the  harsonic  series,  and  any 
finite  sua  of  spherical  harmonics  converges,  indeed  is 
analytic,  everywhere  except  at  the  origin  of  the  coordinate 
systea.  Yet,  if  the  total  infinite  series  of  the  potential 
does  diverge  at  the  earth's  surface,  then  the  aore 
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coefricients  we  determine  the  greater  will  be  the  effect  of 
this  divergent  character  on  the  partial  suas.  The  question 
of  convergence,  however,  is  then  better  posed  as  a  question 
of  representation. 

Since  the  partial  sua,  whether  in  the  space  above  or 
below  the  bounding  sphere,  is  thus  only  an  estimate  of  the 
true  potential,  the  question  of  representation  belongs  to  a 
such  larger  class  of  probleas,  naaely  the  problem  of 
approximation.  As  with  aost  areas  in  physical  geodesy,  our 
liaited  accessibility  to  the  gravity  field,  i.e.  limited  to 
discrete  and  noisy  measurements,  automatically  renders  our 
problem  "ill-posed."  An  ill-posed  problem,  according  to 
Thikonow  and  Arsenin  (1977),  is  a  problem  that  either  has  no 
solution,  has  more  than  one  solution,  or  its  solution  is 
unstable  with  respect  to  the  given  data.  Our  inability  to 
measure  the  gravity  field  in  space  beyond  a  certain  degree 
of  detail  means  that  there  exists  an  infinite  number  of 
solutions,  all  differing  in  the  detail  which  we  were  not 
able  to  discern,  but  all  satisfying  our  measurements. 
Therefore,  the  solution  is  nor  unique.  But  as  easily  as 
this  problem  is  recognized,  it  is  as  quickly  eliminated  by 
requiring  a  solution  for  a  gravity  field  concordant  in 
detail  with  the  measurements .  That  a  solution  always  exists 
is  guaranteed  by  the  fact  that  any  finite  sum  of  spherical 
harmonics,  which  in  fact  satisfy  Laplace*s  equation  and  are 
regular  at  infinity,  represents  a  potential. 

The  instability  of  the  solution  arises  because  the 
harmonic  coefficients  obtained  from  the  measurements  at 
satellite  altitude  are  not  errorless.  This  is  expertly 
shown  by  Bummel  et  ul.  (1979)  for  the  case  that  the 
measurement  noise  is  white  noise.  White  noise  affects  all 
harmonics  of  the  measured  signal  equally  so  that  the 
infinite  sum  of  the  effects  is  unbounded.  Since  the 
discreteness  of  the  measurements  places  a  limit  on  the 
number  of  harmonic  coefficients  that  can  be  determined,  the 
downward  continuation  of  the  error,  while  not  causing 
unbounded  error  in  the  solution,  nevertheless  produces  an 
amplification  of  the  error.  The  error  in  the  n-th  degree 
harmonic  coefficient  is  amplified  by  the  approximate  ratio 
(r/H)n  in  the  process  of  downward  continuation  (see  section 
1.3) ,  where  B  is  the  radius  of  the  earth  and  r  is  the  radius 
of  the  satellite  orbit.  Por  minimal  satellite  altitudes  of 
150  km,  this  ratio  increases  to  over  1000  at  n=300,  which 
means  that  the  300-th  degree  coefficient  of  the  gravity 
anomaly  at  altitude  must  be  known  to  /tgal  (10**  m/s') 
accuracy  in  order  to  recover  mgal  (10*r  m/sv)  accuracy  at 
the  oarth*s  surface.  This  demonstrates  that  the  problem  of 
downward  continuation  belongs  to  the  class  of  ill-posed 
problems. 
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Tae  forekost  method  to  solve  the  ill -posed  problem  is  to 
impose  constraints  on  the  desired  solution,  as  for  example, 
searching  for  a  smoothed  version  of  the  true  solution. 
Another  example  is  the  method  of  collocation  which  generates 
a  solution  whose  norm  in  the  space  of  solutions  is  minimum. 
While  solving  the  nonuniqueness  problem,  collocation  can 
still  be  an  unstable  process,  sometimes  even  requiring  the 
presence  of  noise  in  the  data  to  stabilize,  or  regularize, 
the  solution,  bquivalentl y,  one  can  simply  introduce  a 
regularizing  factor  which  has  the  same  effect  as  noise  in 
that  it  filters  the  higher  frequencies  of  the  solution 
(Ruarnel  et  al.  1979)  .  A  serious  difficulty  with  the  usual 
collocation  is  the  sheer  volume  of  the  computations, 
increasing  with  the  cube  of  the  number  of  the  data. 

However,  with  specially  gridded  data  which  are  then  amenable 
to  very  efficient  computational  algorithms,  Colombo  (1979) 
demonstrates  the  applicability  for  highly  detailed  global 
solutions  of  the  gravity  field. 

Aside  from  collocation,  other  frequently  discussed,  more 
deterministic,  methods  of  downward  continuation  rely  in  one 
way  or  another  on  the  inverse  of  the  solution  to  a 
boundary-value  problem,  either  Poisson's  integral  (first 
boundary -value  problem)  or  the  Pizzett  i-Stok.es  formula 
(third  boundary-value  problem)  ;  both  formulated  on  the 
supposition  of  a  spherical  earth.  When  regarded  as  formulas 
relating  the  sought  after  sources  that  produce  the  observed 
data,  i.e.  as  formulas  for  the  inverse  problem,  they  become 
Fredholm  integral  equations  of  the  first  kind.  Their 
solution  is  usually  found  by  successive  approximations,  but 
because  it  is  unstable,  the  iterations  may  not  converge. 
Assuming  a  spherical  earth,  the  Stokes  "integral  equation" 
is  readily  solved,  yielding  the  inverse  Stokes  equation 
(Holodenskii  et  al.,  1962,  p.50)  .  Host  treatises  on 
applications  of  downward  continuation  were  predicated  on 
airborne  measurements  of  gravity  and  made  use  of  the  Poisson 
integral,  but  only  for  local  determinations;  see  for  example 
the  works  by  Schwarz  (1973)  and  Boritz  (1966a) . 

Finally,  we  note  a  method  of  downward  continuation  that 
is  founded  on  the  usual  technique  for  analytic  continuation. 
Because  the  potential  of  the  earth  (without  atmosphere)  is 
analytic  everywhere  above  its  surface,  the  downward 
continuation  from  the  bounding  sphere  is  theoretically 
achievable  using  a  Taylor  series  expansion.  That  is,  given 
the  potential  on  the  bounding  sphere  (as  a  series)  ,  we  know 
also  its  derivatives.  Hence,  the  Taylor  expansion  about  a 
point  on  the  the  bounding  sphere  can  be  evaluated  anywhere 
within  the  sphere  that  excludes  all  singularities  of  the 
potential;  i.e.  the  sphere  that  is  centered  at  the  expansion 
point  and  just  touches  the  earth's  surface.  This  method  was 
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briefly  developed  by  Uotine  (1969,  pp. 172-173) ,  but  its 
applicability  seems  uncertain.  Continuations  of  3  to  25  km 
using  Taylor  series  require  fairly  accarate  evaluations  of 
tbe  (radial)  derivatives  of  the  potential.  However,  with 
only  a  finite  nuaber  of  haraonic  teras  in  the  series,  the 
derivatives,  particularly  of  higher  ( > 1)  order,  will  suffer 
considerably  froa  the  truncation  effect,  as  well  as  randoa 
errors  in  the  high  degree  coefficients. 


1«3  Preliminaries  and  Definitions 

There  are  certainly  many  additional  aspects  to  the 
problea  of  using  satellite  derived  data  for  terrestrial 
applications.  In  the  first  place,  the  potential  will  not  be 
observed  lirectiy.  The  aeasureaents  at  the  satellite 
altitude  will  consist  of  either  satellite  to  satellite 
Doppler  tracking  data  or  gradioaetry  data.  The  former 
provides  velocity  differences,  hence  potential  differences, 
between  two  satellites  (see  Hajela,  1978;  Buaael,  1980; 
Schwarz,  1972) ,  while  the  latter  yields  linear  coabinations 
of  components  of  the  gravitational  gradient  tensor  (see 
Seed,  1973;  Hummel,  1979).  Secondly,  for  global  coverage, 
the  satellites  must  follow  near  polar  orbits  thereby 
creating  a  nonuniform  data  set  with  heavy  concentrations  at 
the  poles,  in  order  to  perform  a  spherical  harmonic 
analysis,  the  data  must  exist  on  a  sphere  (see  below) ,  but 
the  satellite  orbits  cannot  be  exactly  circular  (the 
satellite  stoves  in  a  noncentral  force  field)  .  Pinally,  the 
earth's  potential  field  is  not  stationary  in  inertial  space. 
It  completes  one  full  rotation  every  24  hours  on  an  axis 
that  wobbles  due  to  precession  and  nutation,  as  well  as 
polar  motion.  Therefore,  the  raw  data  must  undergo 
considerable  preprocessing  in  order  to  obtain  uniform  or 
specially  gridded  coverage  on  a  sphere  that  is  fixed  in  the 
earth's  gravitational  field.  These  preparations  in  the 
determination  of  the  final  product  are  beyond  the  scope  of 
the  following  analysis,  since  they  depend  primarily  on  the 
type  of  satellite  mission,  furthermore,  it  is  assumed  that 
the  potential  is  available  for  downward  continuation  in  the 
form  of  a  (finite)  spherical  harmonic  series.  He  note  that 
the  analysis  of  spherically  distributed  discrete  data  is 
corrupted  by  aliasing,  the  influence  of  the  higher  frequency 
content  of  the  data  on  the  desired  harmonic  coefficients. 
This  effect  can  be  minimized  by  using  optimal  estimation 
techniques  (see  Colombo,  1978)  . 

The  downward  continuation  error  in  the  present  context 
refers  to  a  deterministic  error,  as  opposed  to  a  randoa,  or 
probabalistic  ezror.  Given  the  potential,  defined 
everywhere  in  space,  and  its  series  representation  in  a 


region  of  convergence,  one  can  define  this  error  precisely 
as  the  difference  between  the  truncated  series  continued 
beyond  the  region  of  convergence  and  the  true  value  of  the 
potential.  This  definition  has  the  disadvantage  in  that  it 
includes  a  type  of  truncation,  or  onission,  error,  i.e.  the 
neglect  of  higher  degree  information,  which  has  nothing  to 
do  with  series  divergence.  The  alternative  definition,  as 
the  difference  between  truncated  inner  and  outer  series, 
however,  seems  even  less  agreeable,  since  the  inner  series, 
in  our  case,  is  not  the  spectral  representation  of  the 
potential.  Consequently,  corresponding  terns  of  the  two 
series  are  not  comparable. 

At  present  we  lack  the  resources  (primarily  a 
sufficiently  accurate  series  of  the  potential  determined  in 
outer  space  to  high  degree)  to  conduct  an  analysis  of  the 
downward  continuation  with  actual  data.  The  natural 
alternative  is  to  devise  an  earth  model  with  a  complexity 
that  adequately  takes  into  account  the  anticipated  advances 
in  determining  series  expansions  in  space,  ideally,  the 
potential  of  this  model  should  be  known  exactly  on  the  model 
surface  and  be  expandable  in  a  series  that  diverges  below 
the  bounding  sphere.  Instead  of  exact  values  on  the 
surface,  an  inner  series  may  suffice  if  it  can  be  expanded 
to  a  high  degree. 

From  the  mathematical  standpoint,  the  spherical  polar 
coordinates  r,  9 ,  A  lend  themselves  most  conveniently  to 
formulations  on  a  global  scale.  With  respect  to  the 
Cartesian  system  of  coordinates  x,y,z,  r  is  the  radial 
distance  from  the  origin,  ®  is  the  polar  angle  measured 
from  the  z-axis,  and  *  is  the  angle  (longitnde)  measured 
counterclockwise  in  the  xy-plane  from  the  x-axis: 

x  =  r  sin0  cos  X 

y  ■  r  sin9  sin  X  (1.1) 

z  ■  r  cos0 

in  geodesy,  the  second  coordinate  is  often  the  latitude,  but 
then  is  usually  the  coordinate  in  an  ellipsoidal  system  of 
coordinates.  All  derivations  in  section  2  are  performed  in 
the  above  spherical  coordinate  system. 

The  solutions  to  Laplace's  equation  (which  is  satisfied 
by  the  potential  in  free  space)  are  the  solid  spherical 
harmonic  functions  of  degree  n  and  order  s: 


r“  •  '••<n+1> 


,n>0,-n<m<n  (1.2) 
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the  first  when  the  region  of  harmonicity  contains  the 
origin,  and  the  second  when  it  contains  infinity.  The  Tnm 
are  Known  as  surface  spherical  haraonics  and  are  defined  by 


<e,x) 


fn|m|(cose> 


r  cos  mX 
'sinjraj X 


m  >  0 
m  <  0 


(1.3) 


This  departs  froa  the  aore  conventional  definition  adopted 
by  aatheaaticians  (see  Cushing,  1975,  p.158)  where  the 
inconvenience  of  separate  definitions  for  negative  and 
nonnegative  orders  is  avoided  by  using  the  exponential 
function  e1'"*  instead  ox  the  sinusoids;  however,  (1.3)  is 
nore  customary  in  physical  geodesy.  The  are  the 
associated  Legendre  functions,  noraalized  such  that  the 
integral  of  the  square  of  surface  haraonics_over  the  unit 
sphere  is  4rr.  Furthermore,  the  functions  are 
orthogonal,  i.e. 


1 


//  Ws"x> 
0 


ypq<6.x)do 


r  1  ,  n  =  p  and  m  =  q 
‘0  ,  n^p  or  m#q 


(1.4) 


where  do*  =sin0  d©  dX  and  o'  represents  the  unit  sphere;  and 
they  fora  a  complete  set  of  basis  functions.  This  aeans 
that  any  continuous  function  F(6,X  )  defined  on  the  unit 
sphere,  that  is,  for  0  £  9  £  7T  and  0  £  *  £  2  IT,  can  be 
uniquely  expressed  as  a  uniforaly  convergent  series  of 
spherical  haraonics: 


F( 0 , X ) 


L 

n=0 


n 

I 

m=-n 


nm 


Ynm(0'*> 


(1.5) 


where  (by  multiplying  both  sides  by  YRf<(,  integrating  over 
or,  and  noting  (1.4)) 


t 


nm 


£rff  F<0.*>  ?nm(e *X>dcr 

a 


(1.6) 


The  condition  that  F  be  continuous  can  be  relaxed  to  P  being 
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Lebesque  integrable,  but  then  tbe  series  (1.5)  does  not 
always  converge  to  P.  The  operation  (1.6),  resulting  in  the 
coefficients  f*m,  is  known  as  the  Legendre  transform  and  the 
coefficients  constitute  the  (Legendre)  spectrum  of  P.  in 
view  of  (1.2),  the  extension  of  an  haraonic  function  into 
(the  exterior)  space  is 


F(  r , 0 , X  )  =  l 


I  (|-> 

n=0  m=-n 


n+1 


fnm  Ynm<9>X> 


(1.7) 


where  R  defines  the  radius  of  the  sphere  on  which  P  has  the 
spectrua  {fN»]  (see  equation  (1.6)).  If  the  spectrua  of  P 
is  deterained  in  space,  then  "downward  continuation"  siaply 
means  a  decrease  in  the  variable  r.  Also,  the  spectra  of  P 
on  two  different  spheres,  of  radii  S,  and  are  related  by 


nra  Ri  nm 


(1.8) 


provided  that  the  series  converges  on  each  sphere.  It  is 
obvious  from  (1.6)  that  the  definition  of  Legendre  spectrua 
is  not  restricted  to  functions  defined  on  a  sphere.  The 
surface  can  assuae  any  shape  as  long  as  to  each  coordinate 
pair  (Q,X)  there  corresponds  a  unique  point  of  the 
surface,  and  vice  versa.  Of  course,  if  a  function  is 
defined  in  three  dimensions,  then  its  spectra  with  respect 
to  a  sphere  and  soae  other  (nonspherical)  surface  are  not 
coaparable. 

This  introduction  to  spherical  haraonics  concludes  with 
a  statement  of  a  very  useful  foraula,  the  addition  theorem 
for  Legendre  polynomials: 


P  (cos  Ifl) 
n 


1 

2i.+l 


n 

I 

m=-n 


Ynm<9’X> 


Ynm(e*,X’) 

nm 


(1.9) 


where  P  is  the  central  angle  between  points  (0  ,  A  )  and 
(0',X')  on  the  unit  sphere,  and  where  the  P„  are  the 
faailiar  Legendre  polynomials . 

Vie  follow  Hobson  (1965,  pp. 89-90)  in  the  definition  of 
the  associated  Legendre  functions.  For  any  complex  yu  not 
on  the  real  line  segaent  [-1,lJ, 
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p>>  -  O'1-1’4”  -^T  v> 

«>)  »  (MJ-i)4m  <jn(u) 


(1.10) 


where  the  Q^#  If*  are  the  Legendre  functions  of  the  second 
Jtind.  The  Legendre  functions  with  real  arguaents  M.  =cos  9  £ 
£-1,l]  are  then  defined  by 


Pm(cos  6)  =  (±i)m  p“(cos0±Oi) 
n  n 


(1.11) 

=  (±i)m  lim  P®(cos  e  ±  i/i  ) 
v+0  D 

The  right  side  of  (1.11)  is  the  liait  through  the  coaplex 
plane  onto  tne  line  segaect  Q-l,l]  and  is  found  to  be 


P®(cos0±Oi)  =  (±i)m  sinm0  - - - P  (cos  0)  (1.12) 

n  d(COS0)m  n 


and  siailarly  for  so  that 


p“(cos  0)  =  (-l)10  sinra0  - - -  P  ( cos  0) 

n  d(COS0)m  n 

Q»(cos  0)  «=  (-l)ra  sinm0  - —  Q  (cos0) 

°  d(cos  0)m  n 


Letinitions  (1.10)  and  (1.13)  hold  for  any  n,a,  but  we 
consider  only  those  functions  for  which  n,a  are  nonnegatiwe 
integers  with  0£a—  n.  Finally,  we  apply  the  required 
noraalization: 
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-l)ra 


-1) 


m 


V 


(2n+l)(n-m) ! 
emdn+m) ! 


(2n+l)(n-m) ! 

E_dn+m) ! 

ID 


Q>> 


for  all  y 


for  all  y 


(1.14) 


2.  £h$  pgggfBfd  Con tin nation  2l  Spherical  Harmonic  Series 


Although  the  problea  of  convergence  and  divergence  has 
generally  been  recognized,  lore  attention  has  been  paid  to 
the  theoretical  concerns  than  to  a  noaerical  analysis  of  the 
situation  in  pratice.  1  passing,  conjectural,  coanent  by 
Cook  (19b7)  suggests  that  the  effect  of  divergence  on  the 
n-tb  degree  coefficient  is  on  the  order  of  for  the 
potential  and  J»  for  the  gravity  anoaaly,  vhere  J*  is  the 
n-tn  degree  zodaI  haraoaic  coefficient.  The  notable  works 
in  this  area  are  those  of  Levallois  (1972)  and  Sjoberg 
(1977) .  Their  nuaerical  results  derive  froa  the  postulated 
effect,  oa  the  series  expansion,  of  the  aasses  between  the 
bounding  sphere  and  the  sphere  of  computation,  with  several 
approximations ,  Levallois  estimated  these  effects  over  auch 
of  the  earth *s  surface  for  expansions  of  the  geoid 
undulation  up  to  degree  200.  He  obtained  errors  of  a  few 
tens  of  centime  ter a  with  the  exceptional  aeter  in  eguatorial 
regions,  even  for  low  degree  expansions;  several  meters  in 
the  aidlatitudes,  generally  for  all  degrees  of  truncation; 
and  up  to  16  a  iu  polar  regions  for  the  high  degree 
expansion.  Sjoberg *s  analysis  is  restricted  to  expansions 
up  to  degree  16  and  24,  but  enlarged  to  include  the  errors 
in  gravity  anomalies,  his  results  show  errors  as  large  as 
0.5  to  5  a  depending  on  the  complexity  of  the  earth  model, 
as  well  as  the  point  of  computation .  Downward  continuation 
errors  in  gravity  anoaalies  were  found  to  be  extraordinarily 
large,  in  some  instances  on  the  order  of  the  anoaalies 
themselves  (up  to  30  agal)  . 

These  results  for  the  gravity  anomaly  are  unacceptable 
as  we  have  only  to  coapare  the  GEB9  harmonic  coefficients 
(derived  solely  froa  observations  of  satellite  orbit 
perturbations,  Lerch  et  al.,  1977)  and  coefficients  derived 
froa  terrestrial  data.  Bapp  (1978)  found  an  BBS  (root  aean 
square)  difference  in  the  two  expansions  (up  to  degree  20) 
of  7.0  agal.  He  also  computed  an  BBS  difference  of  9.1  a  in 
the  expansions  of  the  geoid  undulation,  rather  high,  but 
aore  likely  due  to  measurement  errors  than  the  divergence  of 
the  series. 

An  appraisal  of  this  method  of  estimating  the  downward 
continuation  error,  given  in  section  2.2.4,  suggests  that  it 
is  an  unsuitable  aetbod  on  account  of  the  simplistic  density 
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hypothesis  of  the  intervening  Basses.  Also  an  atteapt  is 
made  to  explain  the  irreconcilable  downward  continuation 
errors  of  the  gravity  anoaaly,  aentioned  above.  In  section 
2.3  we  eabark  on  a  siailax  coarse  to  ascertain  the  effects 
of  series  divergence,  bat  now  with  a  firaer  control  on  the 
generation  of  the  earth's  disturbing  potential.  Results 
troa  a  subsequent  nuaerical  analysis  agree  generally  with 
expectations  based  on  the  discussions  of  the  following 
section. 


2.1  Siaple  Hass  Distributions 

A  study  of  potential  series  corresponding  to  siaple  nass 
distributions  will  illuainate  soae  of  the  broader  aspects  of 
the  haraonic  series  behavior  of  the  earth's  potential,  in 
order  to  draw  definite  conclusions  on  the  partial  suas  of 
the  series,  the  Bass  distribution  should  be  sufficiently 
eleaentaxy  so  that  1)  the  surface  separating  regions  of 
convergence  and  divergence  is  well  defined,  2)  the  series 
for  either  region  xs  calculable  to  arbitrarily  high  degree, 
and  3)  although  it  is  not  essential,  a  closed  foraula  of  the 
potential  is  available  to  check  the  nuaerical  coaputations. 

The  following  "experimental  nass  distributions"  are  not 
designed  to  siaulate  the  earth's  distribution  of  Bass; 
hence,  any  of  the  specific  quantitative  results  obviously  do 
not  hold  for  the  earth.  Infinite  series  will  be  developed 
for  both  the  potential  and  the  gravity  anoaaly  since  they 
are  associated  with  different  rates  of  convergence. 

In  the  usual  spherical  coordinate  systea  (1.1),  consider 
an  infinitesimally  thin  layer  of  unifora  density  distributed 
in  the  fora  of  a  circular  disk  on  the  equatorial  plane 
0  s90°  and  centered  at  the  origin  (see  Fig.  2).  The 
constant  density  is  X  and  the  radius  is  denoted  a.  For  any 
point  P,  the  potential  due  to  the  attracting  aass  is 


Vp  =  V(rp,0p,Xp) 


a  211  „  ^  1(.  „ 

t  t  r  d  Ad  r 

<  X  J  1  - 1 - 

r-0  X-0  * 


(2.1) 


where  K  is  the  constant  of  gravitation,  and 

l  »  /rp  +  r2  -  2rrp  cosij>  ;  cos'P  -  sinep  cos(A-xp)  (2.2) 
Because  of  rotational  synaetry,  we  aay  choose  *P*Q.  The 
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Figure  2:  Equatorial  disk  density  distribution. 
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integration  with  respect  to  r  is  easily  performed: 


V  »  k  x  /  [i  -  r  +  r  co&pXn  A.?-ia~rP.co^  ^  ]dA 
p  ■'o  a  p  p  rp-rp  cos  ^  J 


(2.3) 


where 


A&  **  /a*  +  -  2arp  cosifi 


(2.4) 


As  reference,  we  choose  the  potential  of  the  entire  aass 
concentrated  at  the  coordinate  origin: 


U 

P 


*wa2x 

rP 


Hence  the  disturbing  potential  is 


(2.5) 


U 


P 


(2.6) 

<*a2x 


where  because  of  synnetry,  the  integrals  over  the  intervals 
(0,7 T)  and  (77,277)  are  identical. 

The  series  expansions  of  ?  are  obtained  by  substituting 
into  the  integral  (2.1)  the  unifornly  convergent  series  for 
the  reciprocal  distance: 


2k  x  I 


r  +  r 
P  P 


C03»tnY;^f  ]dX  - 


p  p 


oo  n 

*"1  “  l  -j^P^cos*)  ,  r  >  r  (2.7) 

n=0  r“  v 

P 

for  the  outer  series,  and 

JTl  “  l  -^V®08*5  *  rp  <  r  (2*8) 

n™0  r 


for  the  inner  series. 


We  find  for  rf  >  a 


p  “  K  XnL  n+5  (  rD^+2  /**  PnCcos^)dX 


(2.9) 


Now  the  addition  theore.  (1.9)  provides  the  integrals 


of  p„; 


J.2"  '.<“•*>«  -  ^  f  »_«  0)  f2’ 


m»-n  n®'  P'~'  Jo  xnm' 
Tm.  pna  (co SV  Pn/0) 


(2.10) 


since 
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4  Ynm(iir,A)dA  -  0  ,  m  *  0 


Now 


PB<0)  - 


Izl2inl3  5  .  .  .  (n-n 

2®n  (  nj;  i  n  is  even 


n  is  odd 


and  the  series  expansion  for  V  becoies 


(2.11) 


(2.12) 


(2.13) 


where 


un  *  SOT  p2a<°>  i  Pn+1 


The  distorbinq  potential  is  siiplf 


-_E±|_  u 

n+l_  * 


•  n  2  0  ;  w,  *  i  (2.14) 
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2ir 


<X  a  l  (-|“^n+1  u  P,  (cos0, 


Por  points  inside  the  bounding  sphere,  we  decoapose  the 
integral  (2.1)  into  the  two  parts  for  which  the  series  (2.7) 
and  (2.8)  respectively  converge: 


2  <jj  j«_  CO 

vn  =  <x!  t  J  x  (-|-f+1Pn(cosi|Odr  + 
P  n=0  rP  n 


+  /  t  (-^f  Pn(cosi|Odr]dX 
rp  n=0 


(2.16) 


2ir 


“K*rP  J0  /.  tnT3-  PB<co«*)d»  * 

n^l 

2ir 

+  k  X  rp  |9  t-i-  +  tn  —  ]  cos<J>  dX 
The  last  tera  vanishes;  and  with  (2.10)  and  (2.12)  we  obtain 
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The  disturbing  potential,  with  respect  to  the  reference 
potential  (2.5),  is 
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That  this  series  converges  for  rf <  a  is  obvious  once  we 
recognize  that  /^■©(n*1'1)  since  Ptn  (cos  9  )  =0  (n~'/l) . 

Por  the  present  purposes,  we  nay  adopt  the  following 
definition  of  the  gravity  anoaaly  (Heiskanen  and  Boritz, 
1967,  p.89) : 
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Applying  (2.19)  to  (2.6)  and  oaitting  the  tedious 
derivations,  the  closed  fora  tor  the  anoaaly  is  found  to  be 
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Because  surface  layers  generate  potentials  whose  derivatives 
are  discontinuous  as  they  cross  the  surface,  expression 
(2.20)  is  valid  anywhere  except  on  the  disk. 

The  definition  (2.19)  applied  to  the  series  for  the 
disturbing  potential,  (2.15)  and  (2.18),  yields  the  outer 
and  inner  series  for  the  gravity  anoaaly: 
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The  ter»  by  ter*  differentiation  of  the  series  (2.15)  and 
(2.1b)  is  permitted  since  the  series  of  radial  derivatives 
is  uniformly  convergent  with  respect  to  r? . 

Another  siaple  aass  distribution,  which  approaches  that 
of  tne  earta,  is  a  homogeneous  ellipsoid.  In  order  to 
evaluate  the  coefficients  of  the  series  expansions  exactly 
and,  more  importantly,  so  that  the  bounding  sphere  is  the 
surface  of  convergence,  the  ellipsoidal  surface  is  broken 
into  latitudinal  bands  of  constant  curvature,  each  5*  wide, 
giving  it  a  serrated  appearance  (see  Pig.  3).  The  potential 
at  a  point  P  due  to  the  saooth  ellipsoid  is  given  by 
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where  <r  denotes  the  unit  sphere,  da"  =  sin9  d9  dX  ,  and  where 
r,  =is  ( © )  is  the  radial  distance  to  the  ellipsoid  surface, 
with  the  introduction  of  the  serrations. 
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where  ©i=iA0  ,  i=0  ,  .  . .,  2N-1,  H=90/a9  ,  and  rJi=£r1(9^).  Only 
the  inner  integral  can  be  evaluated  in  closed  fora,  thus 
precluding  the  coaputation  of  exact  values  of  the  potential 
and  gravity  anomaly.  The  derivations  of  the  series  are 
coapletely  analogous  to  those  of  the  equatorial  disk  and  are 
relegated  to  Appendix  B.  The  final  results  are  given  by 
equations  (3. 7),  (B.8)  for  the  potential  and  (B.9),(B.10)  for 
the  gravity  anooaly. 

For  the  numerical  tests,  the  equatorial  disk  was  given  a 
radios  of  a=6378l40  a  and  a  uniform  density  of  7c  =3x10*' 
g/cmx.  Similarly,  a =637 8140  a  was  chosen  as  the  equatorial 
radius  of  the  homogeneous  ellipsoid.  The  centers  of  the 
latitudinal  bands,  each  A  9  =5°  wide,  lie  on  an  ellipsoid 
with  a  flattening  of  f= 1/298. 257;  and  therefore,  the 
corresponding  radii  of  these  bands  are  computed  according  to 
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(2.25) 


where  i=1,...,17,  and  e’1,  =1/(1-£)* -1.  The  homogeneous 
density  of  the  ellipsoid  was  equated  with  the  average 
density  of  the  earth,  y.  =5.5  g/ca3 ;  and  K  =66.7x10*’ 
ca/(g.sl)  was  adopted  as  the  constant  of  gravitation. 

The  differences  between  the  resulting  partial  suns  of 
the  series  (2.15),  (2.21),  (B.7J  ,  and  (B.9),  as  functions  of 
the  truncation  degree,  and  the  corresponding  true  values  are 
shown  in  Pig.  4  for  r,  =6377200  a,  9f  =77*5  (near  the 
equator)  and  in  Pig.  5  for  r?=6357200  a,  9,  =7*5  (near  the 
pole).  Both  points  of  evaluation  (rr,9,)  were  selected 
below  the  sphere  of  convergence,  so  that  each  of  the  series 
aust  diverge  in  the  liait.  (Using  Ecker's  theorea  (see 
section  1.1),  it  is  possible  to  prove  that  the  bounding 


sphere  of  the  equatorial  disk,  is  also  the  surface  of 
convergence.  A  siaple  proof  of  the  divergence  inside  the 
bounding  sphere  of  the  serrated  ellipsoid  was  not  found,  but 
it  is  verified  by  Pigures  4  and  5.)  The  partial  suns,  being 
in  any  case  discrete  functions,  were  evaluated  in  steps  of 
20  (Pig.  4)  and  100  (Pig.  5)  degrees  and  connected  by 
straight  line?,  tor  clarity,  but  thereby  also  smoothing  their 
strong  oscillatory  behavior.  The  true  values  (not  shown)  of 
the  potential  and  gravity  anonaly  were  provided  by  fornulas 
(2.6)  and  (2.20)  for  the  equatorial  disk  and  by  the  inner 
series  (B.8)  and  (8.10)  for  the  homogeneous  ellipsoid 
(truncated  at  n=n=30000) .  In  the  figures,  the  value 
immediately  above  each  zero,  related  only  to  it,  indicates 
roughly  the  range  of  the  oscillations  over  the  given  doaain 
of  truncation  degrees. 

A  study  of  these  graphs  reveals  several  interesting 
aspects  of  harmonic  series  divergence.  The  most  obvious 
conclusion  is  that  the  more  distant  the  point  of  evaluation 
is  from  tne  sphere  of  convergence,  the  more  severe  is  the 
divergence  of  the  series.  The  series  near  the  pole  shows 
definite  signs  of  divergence  around  n=1200  to  1800,  while 
the  series  near  the  equator  had  to  be  summed  to  n  >160  00, 
and  higher  for  the  potential,  in  order  to  detect  a 
significant  divergence  pattern.  Also,  for  low  n,  the 
deviations  from  the  true  values  actually  decrease  with 
increasing  truncation  degree  before  they  start  their 
eventual,  unbounded,  increase.  This  is  particularly  the 
case  for  the  potential,  even  in  the  worse  situation  at  the 
pole.  The  distinction  between  truncation  error  and  downward 
continuation  error  (due  to  divergence,  see  section  1.3) 
thereby  becomes  exceedingly  nebulous  for  tbe  lower  degree 
sums.  Whether  a  comparison  of  these  partial  sums  with  those 
of  the  inner  series  gives  a  better  indication  of  the 
downward  continuation  error  is  questionable  for  lower  degree 
expansions,  since  this  type  of  comparison  is  associated  with 
other  interpretive  problems,  as  discussed  in  section  2.2.2. 
The  difference  between  inner  and  outer  partial  sums,  shown 
in  Figures  6  and  7  for  the  sane  two  points  as  above,  does 
indicate  that  some  of  the  truncation  effect  is  common  to 
both.  This  is  particularly  the  case  for  the  equatorial 
disk#  less  so  for  the  ill  behaved  series  of  the  ellipsoid. 
Connected  with  the  overall  delay  in  divergence  is  the 
difference  in  behavior  between  the  potential  and  anonaly 
series.  Since  most  of  the  power  of  the  potential  is 
concentrated  in  the  lower  degree  harmonics,  the  divergent 
character  of  the  series  is  apparent  somewhat  later  than  for 
the  anomaly#  its  power  being  spread  more  to  the  higher 
degree  terms  (i.e.  it  is  essentially  the  derivative  of  the 
potential)  .  These  conclusions  hold  equally  for  both  density 
distributions#  but  are  clearly  more  vividly  depicted  for  the 


Figure  4:  Partial  subs  of  spherical  haraonic  series  of  49, 
T  for  equatorial  disk  and  serrated  ellipsoid 
density  distributions  (Pigures  2  and  3)  sinus 
corresponding  true  values  evaluated  at  r_=6377200 
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equatorial  disk.  The  spikes  in  the  graphs  corresponding  to 
the  hoaogeneous  ellipsoid  are  undoubtedly  due  to  the 
salience  of  the  bounding  surface. 

The  conspicuous  oscillatory  character  of  the  partial 
sums  is  evidently  attributable  to  the  syaaetry  and 
hoaogeneity  of  the  density  distributions.  However,  soae 
type  of  irregular  oscillation  of  a  full  spherical  harmonic 
series  should  not  be  excluded.  Figure  8  shows  the  partial 
suas  of  a  gravity  anoaaly  series  derived  froa  the  (180,180) 
haraonic  coefficient  solution  of  the  earth *s  gravity  field 
(Rapp,  1978)  and  randoa  higher-degree  coefficients  that  were 
scaled  to  decay  according  to  the  Tscherning/Rapp  degree 
variance  aodel  (for  aore  details,  see  section  2.3.2) .  The 
suas  were  evaluated  (using  an  equation  such  as  (2.50))  up  to 
degree  1500  at  a  colatitude  of  10°  (near  the  pole)  for 
various  radial  distances.  Unfortunately,  because  of 
inevitable  constraints  in  coaputer  storage,  and  also  tine, 
such  computations  are  feasible  only  for  polar  latitudes 
where  the  (noraalized)  Legendre  functions  of  high  order  are 
virtually  zero  and  can  be  safely  neglected.  On  the  other 
hand,  since  the  true  value  of  the  anoaaly  below  the 
convergence  sphere  cannot  be  known,  this  graph  is  alaost 
useless  for  quantitative  assessments  of  the  downward 
continuation  error.  The  intent  of  Fig.  8  is  to  illustrate 
the  tendency  for  oscillation  of  a  divergent  series,  as  well 
as  the  very  moderate  effect  of  divergence  tor  degrees  of 
truncation  less  than  300. 


2-2  Th?  Volumetric  Density  Nodel 


2*2.1  The  Derivation  ot  the  Error  Series 

The  masses  of  the  earth  generate  a  potential  wnose 
gradient  is  the  attractive  force  field  postulated  by  Newton. 
It  can  be  shown  that  this  Newtonian  potential  at  a  point  P 
is  the  sua  of  all  attracting  nasses,  each  divided  by  its 
distance  iroa  the  attracted  point  P.  Formulated  for  a 
nondiscrete  nass  distribution,  the  potential  Vp  is 
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K  is  Newton's  gravitational  constant.  Cl  is  the  (three- 
dimensional)  voluae  containing  the  mass  distribution;  is 
tne  density  function;  dH  is  an  element  of  voluae,  so  that 
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the  elemental  mass  is  yw.dQ;  and  X  is  the  distance  between 
the  attracted  point  P  and  the  attracting  nass  /tdfl.  The 
only  restriction  on  the  density  /a  ,  formulated  in  the 
context  of  modern  potential  theory,  is  that  it  be  Lebesgne 
integrable  (Werner,  1974);  hence  discontinuity  of  /a  on  a 
set  ol  measure  zero  is  permitted,  but  the  total  mass  must  be 
finite,  without  significant  loss  of  generality  we  will 
adhere  to  the  classic  requirement  of  piecewise  continuity 
and  boundedness.  The  integral  expression  above  for  V  is 
valid  anywhere  m  space  (i.e.  it  can  be  shown  that  it 
converges  to  the  potential  everywhere,  even  where  the 
integrand  is  singular  (i  =0),  see  (Kellogg,  1953,  p.151)). 
Our  interest  lies  only  on  and  outside  the  surface  that 
bounds  all  generating  masses,  viz.  the  earth's  surface. 

The  expansion  of  V  into  a  series  of  spherical  harmonic 
functions  can  be  founded  directly  on  its  being  a  solution  of 
Laplace's  equation.  Alternatively,  to  give  a  physical 
meaning  to  the  ensuing  coefficients  of  the  series,  the 
potential  \/X  (generated  by  a  point  of  mass  1/k.  )  is  first 
expanded  as  a  spherical  harmonic  series  (equation  (2.7)). 
Upon  the  substitution  ol  (2.7)  into  (2.26),  the  integration 
may  be  performed  term  by  term  to  yield 
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The  validity  of  this  expression  is  guaranteed  only  for 
points  outside  the  bounding  sphere  5h  (see  Fig.  9)  .  In 
(2.2 7)  r  is  the  radius  to  the  volume  element 
dn*rxsin®  d©  dA  dr,  and  cosy  is  now 


cos<p  =  cos9  cos8„  +  sin0  sin6  cos(A-A  ) 
P  P  P 


(2.28) 


The  potential  at  a  surface  point  P  below  the  bounding  sphere 
can  also  be  represented  as  a  convergent  series  by 
considering  separately  those  regions  for  which  the  series 
(2.7)  and  (2.8)  are  respectively  convergent: 
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where  <r  denotes  the  unit  sphere,  dc  *sin6  d6  dX  ,  rs  is  the 
radial  distance  to  the  earth's  surface,  and  r-ain  (rs  ,r,  )  . 
Note  that  rs  is  a  function  of  (  9 , A  ) ,  as  is  r  which  denotes 
the  radius  either  to  the  sphere  of  coaputation  or  to  the 
earth's  surface,  whichever  is  less. 

Using  the  addition  theorem  for  the  Legendre  polynomials 
(1.9),  the  potential  Vp  ,  for  points  above  the  bounding 
sphere,  can  be  expanded  as  follows: 
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where 
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and  where  8  is  the  radius  of  the  bounding  sphere  5).. 
Similarly,  the  expansion  of  the  potential  at  P  inside  Sw  is 
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Bence,  the  expression  (2.32)  is  valid  in  all  of  the  exterior 
space,  he  note  also  that  the  series  (2.32)  tor  rp<H  is  not 
a  series  of  solid  spherical  haraonics,  nor  is  it  an  analytic 
representation  of  V  everywhere  above  the  sphere  of  radius  rp 
because  the  potential  has  discontinuous  second  derivatives 
on  and  inside  the  earth's  surface. 

In  current  practice  the  expression  (2.30)  for  Vp 
(truncated  at  n -h)  is  used  anywhere  above  the  earth's 
surface  (even  inside  Sk) •  The  difference  between  the 
computed  potential  (equation  (2.30)  with  o®  replaced  by  n) 
and  tne  true  potential  (2.32)  is  the  total  error  of 
computation : 
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£.  (?p )  here  is  called  the  downward  continuation  error.  Rost 
authors  identify  only  the  first  tern  as  the  downward 
continuation  error#  in  which  case  the  second  tern, 
representing  the  neglect  of  aore  detailed  infornation,  can 
be  called  the  truncation  error.  However,  for  reasons  to  be 
elucidated  below,  the  first  definition  is  to  be  preferred 
and  will  be  adhered  to  in  all  subsequent  discussions.  Now, 
the  coefficients  of  the  error  series  (2.35)  for  0  &  n  £  n  are 
given  explicitly  by 
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The  first  and  second  terns  differ  by  an  integral  of  the 
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density  in  the  shaded  portion  of  Figure  9;  therefore 


amn<V  =  II  L*  dr 


?  m(6,X,da 
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A  similar  derivation  can  be  found  in  (Sjoberg,  1977;  Cook, 
1967)  . 

In  the  strictest  sense,  the  computation  of  the 
coefficients  a„^(rp),  0  6  n  £  n,  requires  a  knowledge  of  the 
density  of  the  aasses  between  the  coaputation  sphere  sp  and 
the  earth's  surface;  to  assess  the  truncation  error,  we  need 
estimates  of  the  coefficient  functions  (rp)  (or  the 
density  function  of  the  whole  earth) .  f 

Because  the  earth  is  nearly  ellipsoidal  in  shape  and  its 
internal  density,  on  a  large  scale,  exhibits  approximately 
ellipsoidal  symmetry,  the  earth's  gravity  potential,  as  a 
natter  of  convenience,  is  described  with  respect  to  the 
potential  of  a  rotating  equipotential  reference  ellipsoid, 
which  accounts  for  the  coarse  features  of  the  gravity  field. 
How  this  reference  potential,  designated  U,  is  chosen  is 
irrelevant  for  the  problem  at  hand  provided  it  can  be 
calculated  precisely  either  in  closed  fora  or  as  a 
convergent  series  anywhere  on  and  above  the  earth's  surface. 
To  siaplify  subsequent  definitions,  we  also  stipulate  that 
the  potential  on  the  ellipsoid  equals  the  geoidal  potential. 
U  includes  the  aentrifugal  potential  arising  froa  the 
earth's  rotation,  which  therefore  does  not  contribute  to  the 
disturbing  potential.  An  expression  for_U  is  found  in 
(Heiskanen  and  Horitz,  1967,  p.67)  .  If  0  denotes  the 
gravitational  potential  of  the  rotating  ellipsoid  (i.e. 
without  the  explicit  centrifugal  part) ,  then  the  disturbing 
potential  is  given  by 
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since  the  centrifugal  potential  has  been  ositted  in  V.  By 
expanding  the  uoraal  potential  in  a  series,  we  have  from 
(2.30) 
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where  tnm  is  the  difference  between  and  the 
corresponding  coefficients  of  5.  Since  5P  is  supposed  to  be 
known  everywhere,  the  error  in  downward  continuing  a 
truncated  spherical  harmonic  series  of  T  is 


e(Tp)  =  e(Vp)  (2.40) 

where  £(Tp)  is  given  by  (2.35). 

The  aost  ubiquitous  quantity  in  physical  geodesy  is  the 
gravity  anomaly,  simply  because  the  force  of  gravity  is  aost 
readily  observed.  It  is  defined  (when  there  is  no  aass 
external  to  the  geoid)  by 


Ag 


*! 


geoid 


y\ 


ellipsoid 


(2.41) 


where  g  is  the  magnitude  of  the  earth's  gravity  vector  on 
the  geoid  and  y  is  the  magnitude  of  the  gradient  of  0  on 
the  reference  ellipsoid  (see  ibid.  p.  83).  This  definition 
is  easily  generalized  to  gravity  anoaalies  in  the  external 
space  of  the  earth,  where  geoid  and  ellipsoid  are  replaced 
by  geopotential  and  spheropotential  surface,  the  potential 
of  both  surfaces,  in  their  respective  fields,  being 
identical.  Approximating  y  by  the  gravity  of  a  homogeneous 
ball  and  the  normal  gradients  by  radial  derivatives  then 
yields 


(2.42) 


For  a  quantitative  analysis  of  these  approximations,  see 
section  4. 


The  corresponding  spherical  harmonic  series  of  ag  nay 
be  found  as  follows.  Substituting  the  radial  derivative  of 
the  reciprocal  distance  JT*  (see  (2.2)), 


d 

■®p 


j-i  _  -rp+  r  cos  ^ 


into  (2.26)  yields 


(2.43) 


*  W  VTT  (T) 

8  p  * 


I 
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•cJ/J  V 
a 


rp-  r  cos  ^ 


dR 


(2.44) 


Thus,  combining  (2.26),  (2.38),  and  (2.42),  the  gravity 

anosaly  for  any  point  p  becoaes 


Agp-K/;/  yr£i^os!  _  yiL+  a  u 

p  'ft'  43  rpi  3Tp  rp  p 


(2.45) 


Differentiating  (2.7)  with  respect  to  rp  results  in 


^Lrco^Llttl-  .  I  (n+1)  <f  >n+2  Pn(cos*)  ,  r  »  r  ,2.46, 


n=0 


Sisilarly  for  rp <  r,  differentiation  of  (2.8)  with  respect 
to  rp  gives 


rp  ( rp  -  r  cosip)  „  _  j  n(-jp-)n+1  P  (cosij,)  ,  r  <  r  (2.47) 
o*  n=l  r  “  p 


Substituting  the  above  series,  as  well  as  (2.7)  and  (2.8) 
into  (2.45)  : 


dg  -  k  l  [///  y(n-l) Pn(cos^)d«  + 


n=0  (lx 


.n-1 


-  fjj  y(n+2)^g^  Pn(cosi0  d«  ]  -  Gt 


(2.48) 


R2 


where 


Ri  ■  { (r,8,A)  /  rp  >  r}  ,  fl2  *  { (r,6,A)  /  rp  <  r)  (2.49) 

and  where  Gp*-  3  0/  “3  r,  -20/rf  .  »e  note  that  for  points  P 
above  the  bounding  sphere,  the  set  Ol  is  eepty  and  Cl,  s£i; 
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then  with  the  addition  foranla  (1.9)  we  obtain  the  faniliar 
series  expansion 


=  I  l  Hzl  (A)n+2  Y  ( 8_  i  A  )  -  G. 


P  n=0  m=-n  B  rr 


nm  nm  p  p' 


oo  n 


-  l_  l  (f->n+2  g„„  5nm(e„,x„> 


(2.50) 


n=0  m=-n  p 


’nm  nm  p’  p 


by  absorbing  the  expansion  ot  Gp ,  so  that 

n-1 


g 


nm 


R  tnm 


(2.51) 


The  error  in  using  a  truncated  version  of  this  series 
for  points  on  the  earth*s  surface  (i.e.  below  the  bounding 
sphere)  is  therefore  the  difference  between  equation  (2.50) 
(oo  replaced  by  n)  and  the  true  expansion  (2.48): 


rs 


e(Ag.)  -  k  t  JS  r(n-l)  /_~M(p-)n+2dr 

p  n=0  o  r  p 


+  (n+2)  /_S  y(^£.)n-1dr]  Pn(cos^)  da 
r 


-  k  l  SI  [(n-1)  /  u(i) 

n=fi+l  a  #  p 


r  .n+2 


dr 


(2.52) 


-  (n+2)  j_S  p(^)n~1  dr]  P  (cosi|>)  da 
r  r  n 


Then  with  the  addition  theoren, 

n  n 


e(Ag  )  “II  dnm(rp^  *nra(0p,Xp^ 


n=0  m**-n 


♦  I  I  <Wrp>  WVV 

n-fi+i  m«-n 


(2.53) 


f 


I •»  •* 
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where 


d  (r) 
nm  p 


and 


2n+T 


//  J_S  [(n-l)(^-)n+2  +  (n+2)(^.)n-1] 


a  t 


ydr 


(2.54) 


do 


,R  Nn+2  n-1  .  , 

(“)  ~»s  1  v  +  d  (r  ) 

rp  R  nm  nnr  p' 


(2.55) 


(2.55)  can  also  be  verfied  by  adding  and  subtracting  the 
density  integral  over  the  region  .  We  note  that  for 
rp  >  B,  d^tr,,)  =0  and  the  renaining  error  is  simply  one  of 
truncation. 


Similarly,  the  geoid  undulation  (or  aore  generally,  the 
height  anomaly,  i.e.  the  separation  between  geopotential  and 
spheropotential  surfaces  at  the  same  potential)  is  given  by 
Bruns*  formula  (Ueiskanen  and  horitz,  1967,  p.85): 


(2.56) 


where  y,  is  the  normal  gravity  at  the  point  Q,  being  the 
normal  projection  of  P  onto  the  spheropotential  surface. 

Yq  is  conventionally  approximated  by  its  average  value  on 
the  reference  ellipsoid.  Here,  we  use  a  common  alternative, 
namely  the  gravity,  at  P,  associated  with  an  homogeneous 
ball:  mKH/Vf  (M  *  the  total  mass  of  the  earth)  ;  see 

section  4.  We  then  have 


c  ■ 

^p  kM 


(2.57) 


and  for  points  P  above  the  bounding  sphere 


-  I  1  <£ > 

n=»0  m»-n  P 


R  xn-l 


z  Y  (  9  ,  X) 
nm  nm  p  p 


(2.58) 
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where 


Jnm  k  M  nm 


(2.59) 


Using  only  the  terms  up  to  degree  n  to  represent  the  height 
anomaly  at  the  earth's  surface  results  in  an  error  given  by 


e<Cp)  - 


-7m  e(V 


(2.60) 


where  £(Tp)  is  given  by  (2.40) 


2.2.2  The  Interpretation  of  the  Error  Series 

The  proper  understanding  of  any  nnaerical  computations 
of  errors  such  as  (2.35)  or  (2.53)  cones  only  with  the 
correct  interpretation  of  the  true  series  expansion,  such  as 
(2.32)  or  (2.48).  Consider,  for  example,  the  potential 
series  (2.32) ;  the  sane  arguments  obviously  hold  for  the 
gravity  anomaly.  Recalling  that  B  is  the  radius  of  the 
bounding  sphere  Sb,  the  potential  in  the  space  exterior  to 
Sfc  is  given  by  the  uniforaly  convergent  series  (2.30).  The 
coefficients  vnwt  h*ve  either  of  the  following 
interpretations.  First,  they  are  density  integrals,  as 
given  by  (2.31) .  Secondly,  they  constitute  the  spectrum  of 
V  on  the  bounding  sphere  Sb  (cf.  (1.6)): 


vnn.  “  K  //  VCR, 9 , X)  ?nn,Ce,X)do  ,2.61) 

a 


The  expansion  of  the  potential  in  spherical  harmonics  at 
points  below  the  bounding  sphere  is  achieved  by  considering 
separately  the  two  domains  in  which  the  inner  and  outer 
series  of  solid  harmonics  converge.  The  resulting 
"coefficients"  v^ (rp)  are  more  correctly  functions  of  r, 
and  in  the  first  place  are  density  integrals  (eguation 
(2.33)).  The  vHw,  do  not  describe  the  spectrum  of  the 
potential  on  the  earth's  surface.  That  is,  the  potential 
spectrum  of  constant  coefficients  obtained  from  surface 
values,  as  suggested  by  (1.6),  and  the  functions  vHM>(rp)  are 
clearly  not  identical.  Instead,  because  the  series  (2.32) 
converges  everywhere  to  the  potential  V,  even  inside  the 
earth's  body  (Kellogg,  1953,  p.151),  the  functions  v„m(Bp) 
for  constant  rp=Bf  represent  the  spectrum  of  the  potential 
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on  the  sphere  of  radios  H f ,  whether  inside  or  outside  the 
eaith*s  body;  that  is. 


WV  =  K  i'/v<Rp'9-X>  *>d°  (2-62) 

a 


The  tera  by  tern  evaluation  of  the  downward  continuation 
error,  as  given  by  equation  (2.37)  ,  therefore  does  not 
produce  negative  corrections  to  the  spectrua  of  the 
potential  on  the  surface.  Me  are  forced  to  identify  the 
downward  continuation  error  with  the  entire  sub  of  error 
coaponents  and  on  a  point  by  point  basis. 

Obviously,  evaluations  of  the  infinite  sub  are  beyond 
our  coaputational  ability.  Moreover,  the  finiteness  of  the 
nuaber  of  haraonic  coefficients  of  the  potential  deterained 
in  space  necessarily  liaits  our  efforts  to  estiaating  Bean 
values  or  the  potential  or  gravity  anoaaly,  and  not  point 
values.  This  requires  a  aodification  in  the  foraulation  of 
the  downward  continuation  error.  Consider  first  functions 
defined  on  a  sphere  and  define  an  (isotropic)  averaging 
operator  by 


F(  0  ,  X  )  =  //  B ( \p )  F(  0  ,  X  )do  (2.63) 

% 

where  P  denotes  the  average  of  a  function  P  over  the 
circular  cap  <rc=  {(fl,A)/0  SY  ‘T(]  and  weighted  by  the 
kernel  B(y').  y-is  the  angle  between  the  center  of  the  cap 
(©,>)  and  the  point  of  the  integration  (0,A).  fteissl 
(1471,  p.26)  sh ows  that  if  the  kernel's  (one  diaensional) 
spectrua  is  {  \Ji n*  1  /3„/4Tr]  and  the  spectrua  of  P  is  {f.J, 
then 

fnm  s  Bn  fnm  (2.64) 


are  the  spectral  coaponents  of  the  average  of  P,  that  is,  P. 
In  the  tera  os  o  logy  of  spectral  theory,  the  coefficients  p„ 
are  also  known  as  the  frequency  response  of  the  averaging 
operator  (2.63) . 

Applying  this  result  in  the  present  context,  let  7 
denote  a  truncated  version  of  the  general  series  (2.32) , 
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valid  anywhere  on  or  outside  the  earth's  surface: 


.  n  n 

Vp  -  l  l 

H  n=0  m=-n 


*nm<rp>  Ynm<VV 


(2.65) 


Por  constant  rp-Bp,  this  nay  be  interpreted  as  a  weighted 
average  ot  v  over  the  sphere  of  radius  Bp ,  where  the 
frequency  response  is  unity  for  0  4  n  4  n  and  zero  for  higher 
degrees.  It  is  impossible  to  devise  such  a  perfect  response 
for  the  usual  average  that  is  Halted  to  values  within  a 
cap.  Ve  consider  instead  the  average 


<X> 


n 


vnm(rp> 


Ynm<VV 


(2.66) 


where  the  weights  have  been  chosen  so  that  the  freguency 
response  is 


a  0  ,  n  >  n 
n 


(2.67) 


a  value  to  which  it  tapers  suoothly  frou  a  value  of  1  at 
n=0.  is  an  average  on  the  sphere  of  radius  tr ,  it  is  clear 
that  values  of  V  coincide  with  values  of  the  potential 
averaged  over  an  area  ot  the  earth's  surface  only  if  that 
area  coincides  with  the  spherical  cap.  This  is  never  the 
case  exactly,  but  it  is  a  reasonable  approximation  if  the 
cap  is  swall.  The  size  of  the  cap  for  the  average  (2.66)  is 
essentially  deterained  by  the  desired  "cut-off  frequency"  n. 
Therefore,  if  the  potential  as  deterained  outside  the 
boundinq  sphere  is  first  averaged  over  a  spherical  cap,  or 
equivalently,  its  spectrum  is  aultiplied  by  a  particular 
frequency  response  function,  and  the  truncated  series  is 
evaluated  at  the  earth's  surface: 


A 


n 

n-0 


n  n 

I  <f> 

m*-n  p 


n+1 


Jn  vnm  ?nn,(VV 


(2.68) 


then  the  downward  continuation  error  of  the  average 
potential  ?  is  (cf .  (2.35)  and  (2.36)) 
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£(Vp> 


n 

n=0 


n 

l 

m=-n 


anm<rp> 


Ymn<V 


V 


(2.69) 


k  similar  result  holds  for  lean  gravity  anomalies,  as  well 
as  aeao  height  anonalies: 


fCAgp) 


n 

I 


n 

I 


n=0  m=-n 


fi_  d(r)  Y  (e  ,X  ) 
n  nm  p  nm  p  p 


(2.70) 


e(Cp> 


2  n  n 

— Rr  1  1  6.  a__(r  )  Y(0,X_) 

n£0  m=-n  n  nm  P  nm  p  p 


(2.71) 


where  d„„,(r*)  is  given  by  (2.54)  and  in  each  case  the 
frequency  response  /3„  is  assuned  to  vanish  for  n>n 
(equation  (2.67)). 

An  operator  which  filters  higher  degree  harnonics  nearly 
perfectly  is  the  Gaussian  filter,  its  nane  deriving  froa  the 
shape  of  the  weighting  function  wfi  (>*✓),  defined  by 


Bg(UO 


wg4> 

//wG4)do 

a 


(2.72) 


where 


WQ(l|))  = 
s 


e-a(  1-cos  ip) 
e-^8^  ,  for  small 


'P 


(2.73) 


Zf  we  define  w^(^>)s0  for  v  >  ve  »  then  the  frequency 
response  is  given  by  (see  Jekeli,  1981) 


Go 


=  1 


Gi 


l-y.e-a<1-y,> 


i 

a 
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FI 


2n+l 


“Gn+1 


+ 


e-a(l-y9) 
n-1  1_e-a( i-yo 


(2.74) 


* lPn-l<y°  y  "  Pn+l(y®  >]  »  n  >  0 


where  y0 =cos  y; .  The  parameter  "a"  specifies  the  amount  of 
smoothing.  If  we  desire  fi*=f0  (e.g.  f„=0.05  implies  that 
only  5%  of  the  S-th  degree  harmonic  coefficient  contributes 
to  the  average) ,  then  an  approximate  formula  for  "a"  is 


(2.75) 


2.2.3  The 


of  the  Earth  Hod el 


To  evaluate  the  density  integrals  of  the  downward 
continuation  error  requires  a  simplification  of  the  earth's 
surface,  as  well  as  the  volumetric  density.  In  the 
simplifying  scheme  adopted  here,  the  volume  between  the 
sphere  of  computation  and  the  earth's  surface  is  partitioned 
into  (three-dimensional)  blocks  that  are  delimited  on  the 
sides  by  6 -constant,  1  =constant  and  on  the  top  and  bottom 
faces  by  r(=oonstant,  rF  “constants  Hi  thin  each  block  the 
density  is  assumed  to  vary  only  as  a  step  function  in  the 
radial  direction  and  be  otherwise  constant.  The  required 
lithospheric  (crustal  and  upper  mantle)  densities,  as  well 
as  the  corresponding  depths  are  taken  from  (Bomford,  1971, 
p.457)  and  illustrated  in  Pig.  10.  With  these  assumptions, 
the  coefficients  dnm(r,)  of  the  downward  continuation  error 
of  the  gravity  anomaly  (2.54)  reduce  to 


dnm(rp)  *  i»p(  2nn 


l  WVV  II  <2-16' 


where 


Gn(rp,9i’Xi) 


/  Si  [<n-l)(-L)n  +  <n+2)<£B)n+1]ur2dr  , 
r«  rP  r 


rsi  >  rp 
0  •  rsii  rp 


(2.77) 
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Pigote  10:  The  voluaetric  density  aodel  (density  values  ace 
in  units  of  g/cm7)  . 
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Defining  for  rfc  >  ra 

F  (r  ,r.  )  =  f  b  [(n-l)(^-)D  +  (n+2)(^?)n  1]r2dr  >  n  i>  0 
nab  J  ”p  1 


(2.78) 


Sr -  <^>n+3’ '  Si  ^<£>"~2-  <£>"'2i  • 


P  rb  *a 

n  >  0  ,  n  ^  2 


4-  r’t(^)5  -  (Ja)5  ]  +  4r*  £n  ^  ,  n  = 

o  P  rp  rp  p  ra 


the  functions  G„  are  given  by  the  following  linear 
combinations  (see  Fig.  10) 
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►  *\> 

+  W3Fn(r1  ,i*)  +  PiFn 
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o, 

rp<r2: 
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> 

o. 
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p  si 

V 

P>Fn 

<v 

rsi> 

(2.79) 
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0, 

rp<r,: 

Gn= 

“'W 

.rs) 

+  P*Fn 

(r»,r*)  +  y2Pn 

(rH, 

rsi) 
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< 

0, 

rj  < 

rp<r-: 
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P»Fn 
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< 
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<v 
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If  p  denotes  the  geocentric  distance  to  the  sarface  of  a 
reference  ellipsoid,  then  the  radii  are  computed  as  follows 


+  +  ®i  '  >  0 

Pj_  *  i  0 


■****■  -  i»  hM»  ,  ' 
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Ti 

rsi 

— 

Hi 

-  - 

r2 

=  ri 

+ 

20 

km 

r3 

rsi 

- 

11 

km 

Tk 

rsi 

+ 

H. 

l 

where  H;  is  the  geoid  undulation  and 

a 


✓  1  +  e cos^B. 


H1  ^  0 
Hi  <  0 


(2.80) 


(2.81) 


a  is  the  seaiaajor  axis  of  the  ellipsoid  and  e*  is  the 
second  eccentricity,  related  to  the  flattening  f  by 
e,r  =  1/(1-f)^ -1 .  The  surface  of  the  earth  was  divided  into 
latitudinal  bands  each  5°  wide,  and  each  band  was  further 
subdivided  into  blocks  having  approxinately  the  area  of  a 
5°x5*  block  at  the  equator  (Hajela,  1975).  The  elevation 
data  that  Sjoberg  (1977)  used  provided  values  of  H;  (on  land 
or  sea,  see  Pig.  10)  at  the  center  points  (  6i  ,  Xi  )  of  these 
blocks.  Finally,  the  undulation  H,  was  coaputed  fron  the 
GEB10B  haraonic  coefficients  (complete  to  degree  and  order 
36,  Lerch  et  al.,  1978): 


(2.82) 


where  R  is  the  Beau  radius  of  the  earth  (R=6371  ka) . 

The  aodeling  of  the  density  as  above  does  not  confora  to 
any  established  theory  of  isostasy.  Since  graviaetric 
evidence  indicates  that  the  nass  excesses  and  deficiencies 
near  the  earth's  surface  are  to  soae  extent  isostatically 
compensated  deeper  within,  a  model  which  incorporates  this 
idea  aay  be  more  authentic.  Adopting  the  iiry-Heiskanen 
theory  of  isostasy  (Heiskanen  and  Horitz,  1967,  pp. 135-6) 
with  an  assumed  crustal  thickness  of  0*30  ka  changes  the 
density  aodel  only  in  the  way  the  radii  r,  ,  ...,r,,  are 
computed  (changing  also  /l,,  to  3.27  g/ca’ )  : 


ri  -  r2 
r  j  - 
r«,  * 


—  D  —  5 .45  , 

rsi  “  D  “  2.73  Hi j 


>  0 

H.  <  0 


(2.83) 
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2-2.4  The  Nuaerical  Analysis 

He  consider  only  evaluations  of  the  downward 
continuation  error  in  gravity  anoaalies  because  they  support 
aost  vividly  the  subsequent  conclusions.  The  dependence  of 
the  error  (2.53)  on  the  radial  distance  r,  necessitates  a 
point  by  point  evaluation.  A  dense  grid  of  coaputed  values 
over  the  entire  earth  is  prohibited  by  the  excessive 
coaputer  tiae  that  would  be  required.  Therefore,  in  order 
to  estiaate  the  downward  continuation  error,  we  aust 
restrict  ourselves  to  a  judicious  selection  of  points,  which 
should  be  governed  by  our  objective  to  detect  the  influence, 
not  only  of  the  earth's  ellipticity,  but  also  its 
topography.  For  exaaple,  consider  a  single  profile  in 
longitude  across  southern  Africa  (  0f  =1Q2?5,  -2^5  i  A,<  7?!5)  . 
Figure  11  shows  the  earth's  shape  in  this  profile  as  defined 
by  the  5°  mean  elevations  and  the  GEM10B  geoid.  For  each  of 
six  points  along  the  profile  the  first  tera  in  (2.53)  was 
calculated  with  n=16,  36,  and  180.  Fig.  11  shows  that,  as  a 
supposed  downward  continuation  error,  its  values  for  low 
truncation  degrees  are  inordinately,  in  fact  unbelievably, 
large  (see  the  introduction  to  chapter  2) .  As  n  increases, 
this  "error"  generally  decreases  but  not  aonotonically  as 
shown  in  the  Fig.  12  of  partial  suas.  The  sunaation  to 
degree  180  is  not  strictly  legitimate  for  an  earth  that  is 
saapled  on  a  5°  equal  area  grid  of  a  total  of  1654  values. 
That  is,  1654  bits  of  information  on  the  gravity  field 
deteraine  a  aaxiaua  nuaber  of  1654  s  (39+1>x  coefficients  in 
its  spectral  haraonic  representation  (see  also  Shebalin, 
1980) .  Therefore,  the  coaputed  teras  for  degrees  greater 
than  39  in  no  way  reflect  the  earth's  true  gravity  field, 
but  the  error  teras  to  degree  180  are  included  to  illustrate 
their  general  trend.  Of  course,  in  the  evaluation  of  the 
error  according  to  equation  (2.53) ,  we  have  totally 
neglected  the  higher  degree  contribution  froa  5+1  to 
infinity  because  it  is  unknown  for  this  aodel. 

Froa  the  few  nuaerical  results  presented  in  Figures  11 
and  12  and  on  the  basis  of  the  experiaents  described  in 
section  2.1,  as  well  as  the  discussion  in  section  2.2.2,  the 
following  inescapable  conclusion  is  asserted.  The  errors 
depicted  here,  instead  of  showing  the  divergent  character  of 
the  series,  nay  rather  be  a  reflection  of  the  iaplicit 
choice  of  aodel  for  the  earth's  gravity  field  at  the  earth's 
surface.  In  the  first  place.  Figures  4  and  5  of  section  2.1 
strongly  suggest  that  the  divergence  is  not  aanifest  for  low 
values  of  the  truncation  degree,  especially  near  the 
equator.  If  we  accept  this,  then  the  first  part  of  equation 
(2.53),  i.e.  that  being  evaluated  here,  aust  be  the 
difference  between  two  entirely  incoapatible  partial  suns. 

In  fact,  the  downward  continued  sub  represents  what  could  be 


Figure  11:  Topographical  profile  of  southern  Africa 

(9,*  10215,  -12*5  A  90*0)  and  values  of  partial 
suss  of  downward  con  tin  a  at  ion  error  in  gravity 
anosaly  (voluaetric  density  aodel) :  5=15,36,180. 


ial  suas  of  dounvard  continuation  error 
(volumetric  density  «odel)  for  pouts  1 
iuure  11  and  Of  n  £.  100. 
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called  a  "free  air"  spectra*  of  the  gravity  anoaaly,  while 
the  "true*  partial  so*  represents  the  spectra*  of  a  field  on 
a  sphere  that  is  partially  eabedded  inside  the  earth.  These 
arguments  and  the  supporting  nuaerical  analysis  reinforce 
the  statement  Bade  in  section  2.2.2  that  the  individual 
terns  of  equation  (2.53) ,  or  in  the  strictest  sense*  any 
finte  aggregate  of  terns,  do  not  yield  corrections  to 
corresponding  downward  continued  tens .  This  is 
reeaphasized  here  because  just  such  a  procedure  is 
occasionally  iaplied  in  the  literature  (Cook.,  1967; 

Morrison,  1970) . 

It  nay  be  noted  that  the  kernel  of A the  error  integral 
(2.54)  can  be  expanded  as  a  series  in  h=1-r/rp .  They  the 
infinite  series  of  the  constant  and  linear  terns  in  h  of  the 
downward  continuation  error  sron  to  zero  (see  section  2.3.3) . 
Without  these  terns  the  nuaerical  values  of  Pigure  11  would 
decrease  by  1,  2,  or  nore  orders  of  nagnitude.  However, 
this  does  not  alter  the  essential  conclusion  drawn  above. 

If  n  oo  ,  then  the  "true"  partial  sun  converges  to  the 
actual  surface  value  of  the  gravity  anoaaly,  approaching  the 
downward  continued  series  before  it  diverges.  Sane 
indication  of  this  is  given  by  Pig.  12  which  shows  an 
overall  decrease  in  nagnitude  of  the  error  with  increasing 
n.  Thus,  the  next  step  in  the  analysis  would  be  a 
densif ication  of  the  grid  on  which  the  elevations  and 
densities  are  assigned,  thereby  allowing  expansions  of  the 
error  to  higher  degrees.  However,  the  nodeling  of  the 
disturbing  potential  to  a  high  degree  by  voluaetric  density 
distributions,  or  (what  is  alnost  equivalent)  point  aasses, 
on  a  global  scale  is  generally  associated  with  a 
considerable  conputational  effort  (Needhaa,  1970),  as  a 
distribution  of  nasses  is  sought  vhich  fits,  in  a  least 
squares  sense,  our  knowledge  of  the  exterior  gravity  field 
(see  also  Balnino,  1972).  In  our  case,  no  infornation  on 
the  earth's  gravity  field,  except  postulated  nean  densities, 
was  used  to  define  the  distribution.  One  should  therefore 
not  expect  this  type  of  nodel  to  produce  a  close  resenblance 
of  the  earth's  potential.  Even  the  use  of  the  isostatic 
nodel  of  airy  and  Heiskanen  (equation  (2.82)),  instead  of 
the  nodel  depicted  in  Pigure  10,  does  not  produce 
significant  changes  in  the  nuaerical  results.  This  nodel  is 
therefore  not  further  used  in  the  analysis  of  the  downward 
continuation  error.  In  the  following  section  the  earth 
nodel  is  also  deterained  by  a  density  distribution,  which 
although  not  optinal,  is  entirely  adequate  to  produce  a 
realistic  potential. 
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2.3  The  Density  Layer  Ho del 


2.3.1  The  Deriyation  of  the  Error  Series 

Other  than  voluaetric  Basses,  Newtonian  potentials  are 
generated  as  well  by  sin pie  layer  and  doable  layer  density 
distrioutions.  The  fornalation  of  the  potential  of  a  sinple 
layer  is  coapletely  analog  on  s  to  (2.26),  except  that  the 
integration  is  over  a  surface  instead  of  a  voluae,  and,  of 
course,  the  nass  is  distributed  as  an  infinitesiaally  thin 
layer  on  a  surface.  This  leads  directly  to  the  Molodenskii 
boundary-value  problen  where  the  density  of  the  layer,  as  an 
unknown  guantity,  is  related  to  the  boundary  values  of  the 
resulting  gravity  field  through  an  integral  equation.  This 
integral  equation  can  be  solved  readily  if  the  surface  is  a 
sphere  and  by  successive  approxiaations  for  lore  coup  Heated 
surfaces  such  as  the  earth's  surface.  Por  the  present 
purposes,  the  choice  of  foraulation  of  the  disturbing 
potential  is  dictated  by  our  objective  not  only  to  find  a 
reasonable  solution  to  the  density,  but  also  to  expand  the 
potential  in  spherical  haraonic  series  above  and  below  the 
bounding  sphere.  The  solution  for  the  density  will  be 
deterained  appr oxiaately  by  our  knowledge  of  the  gravity 
field,  that  is,  the  bounding  values;  hence  the  disturbing 
potential  generated  by  this  density  layer  should  aore 
faithfully  represent  the  actual  exterior  potential  of  the 
earth.  Note  that  in  the  following,  no  atteapt  is  aade  to 
find  the  optiaal  density  that  fits  our  knowledge  of  the 
gravity  field. 

This  aethod  of  foraulating  the  downward  continuation 
error  originates  with  Petrovskaya  (1979)  .  It  is  here 
derived  in  aore  detail  and  froa  a  slightly  different  angle 
of  approach. 

Let  us  then  consider  the  following  foraulation  of  the 
disturbing  potential  (with  respect  to  a  suitable  reference 
potential),  according  to  Brovar  (1964), 


v(0,X)  E  do 


(2.84) 


where,  as  before,  B  is  the  radius  of  the  bounding  sphere;  v 
is  a  generalized  density  which  contains  also  the  inclination 
tera  that  transforas  the  integration  over  the  earth's 
surface  to  an  integration  over  the  unit  sphere;  and  the 
kernel  E  is  a  function  of  0,  A,  »>,  ®p,  XP,  and  is  defined 


» 


I 
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by 


where  rs=rs  (0,A)  describes  the  radial  distance  to  the 
earth's  surface.  In  general,  the  definiton  of  the  Kernel  is 
contingent  only  on  the  requirement  that  the  function  T  be  a 
potential,  i.e .  harmonic  in  free  space  and  regular  at 
infinity,  otherwise  it  is  arbitrary.  In  order  for  T  to  be 
harnonic,  E  as  a  function  of  the  point  P  Bust  satisfy 
Laplace's  equation  and  be  regular  at  infinity.  This  was 
shown  by  Brovar  (1964)  to  be  the  case  for  the  above 
definition;  in  tact,  the  difference  between  E  and  the 
generalized  Stokes  function  S  is  easily  recognized  to  be 
(Heiskanen  and  Aoritz,  196  7,  p.93) 


qy  o  —  5rscos» 
3E'S~  rj 


(2.86) 


Stokes'  function  and  the  terns  on  the  right  side  of  (2.86) 
all  are  harnonic  functions.  Unlike  the  usnal  density 
integral  in  which  \/l  is  the  kernel,  the  expression  (2.84) 
has  continuous  derivatives  when  crossing  the  surface  and  is 
therefore  characteristic  of  the  potential  of  a  voluaetric 
density  distribution  (see  Brovar,  1964,  and  Boritz,  1966b, 
p.55) . 

The  radial  derivative  of  T  is 


4r-“  *5?//  v<e'x>  da  <2-87> 

pa  p 

where 

3E  -  2j-  rs  co_s  \jj  1  2rscos4>  ,  rp  +  A  -  rs  cos<{>  (2.88) 

8rp  r*  f,rp  rj  2rp 

This  conbined  with  2h/rp  results  in  the  sisple  expression 
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rs  COSti  . 

- “2 -  J 


(2.89) 


For  points  P  on  the  earth's  surface,  we  have  the  boundary 
condition  (with  the  spherical  approximation  as  described  for 
(2.42)) 


(2.90) 


so  that,  because  the  derivatives  of  the  right  side  of  (2.84) 
are  continuous  at  the  surface. 


AgP  "  5^  //  v(6  ,1  >(y  -  go  (2.91) 


and  the  reason  for  the  choice  of  E  is  now  apparent.  The 
expression  (2.91)  is  a  Fredhola  integral  equation  of  the 
first  kind  and  if  the  earth*s  surfaceis  approximated  by  a 
sphere  with  radius  &  (r4  *const*B,  r,  =R)  on  which  aq  is  a 
known  function,  the  solution  for  the  density  is  found  by 
expanding  1/1  and  ag  in  series  of  spherical  harnonics 
(Petrovskaya,  1979;  see  also  equation  (2.7)): 


nl0  m--n  8nm  ?nm(  W  =  Agp 


*  «  n  .  (2*92) 

4ibj0  J  n  !J  V<9,U  ?n*<9’X>do  *mn<VV 

n=0  m=-n  o 

n*l 


This  holds  for  arbitrary  points  P  on  the  sphere,  hence 


1 _ R* 

Sn+T  4¥F 


!!  v(6  ,X  )  ?nm(9,X)do 


(2.93) 


Therefore,  the  spectrua^of  the  density  in  spherical 
approxiaation  is  (2n*1)g„m: 


58 


d  “  n 

v(e,x)  =  (-£■?  l  l 


n=0  m*-n 


<2n+1>  Km  W9'1’ 


(2.91) 


Because  the  gravity  anonaly  has  no  first  degree  ter*  in  its 
series  representation  (equation  (2 .51)),  this  ter*  for  the 
density,  hence  for  the  potential  T,  *ust  be  obtained  fro* 
other  data. 

The  series  expansion  of  ag  follows  directly  fro*  (2.91) 
if  the  boundary  condition  (2.90)  is  used  as  the  definition 
of  the  gravity  anoaaly  in  all  of  space.  Using  (2.7)  and 
(2.8),  we  obtain  for  r,  >  fi 


Ag_  -  l  l  (^-)D+2  gnm  ?  <6  .A) 

p  a±0  m--n  rp  nm  nm  P  P 


(2.95) 


wh  ere 


g, 


nm 


»  n  «  1 


1  1 


l  SHTI  TjlJ  v(9,A)(^f)n  ^nnCe.AJda  ,  n  *  1 


(2.96) 


and  for  rp<  B 


os  n 


AgP  ■  Jo  J-n  “-‘V  V  W  <2-”> 


where 


cnm(r_) 
nm  p 


(  T  5?r3  //  v(9'X>  3  *lro(e»x)do 

P  o2  s  P 


□  «*  1 


’drrwp  l!S  ?nm<9,xxip  +  <2.98) 

’  XP  Oi  P 

^  +  //  v<9,A)(jjj)n+1  ^nmC6,A)do]  ,  n  *  1 
a  2  8 
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•here  <rx  ~  {  (  ©  ,  X  )/r,  (  9  ,  X  )  (.  rt>  and 

<rx  -  { { 0  »  X  )/r,  (  ©  ,  X  )  >  r,}  .  The  downward  continuation 
error  is  the  difference  between  the  truncated  series  (2.95) 
evaluated  at  the  earth's  surface  and  the  true  series 
expansion  (2.97): 


n  n 

‘  Jo  J-n  d-<rP)  VV  + 

+  jL  in  d-<rP)  ’-‘V  V 


(2.99) 


where 


d  (r  )  -  (^-)n+2  gnm  -  c  (r  ) 
nm  p  rp  "nm  nm  p 


•jsrsnr//  ',(9’X)  <2-10<» 


,?nm(e»X)da  *  n  ~  0 


and 


d’  (r  ) 
nrav  p' 


~c(r) 
nm  p 


nm 


CO  -  (??-)n+2 


e 


nm 


(2.101) 


The  expansion  into  series  of  the  disturbing  potential  is 
sore  involved  and  requires  the  expansion  of  the  kernel  E 
into  series  which  are  valid  below,  as  well  as  above,  the 
bounding  sphere.  Anticipating  the  result,  let  us  consider 
the  integral  (tiradshteyn  and  Byzshik,  1980,  p.83) 


/  drp  -  l  +  r#  cost/;  In  2<rp+t  -  rg  cos<M  +  a(r*)  (2.102) 


where  a(r$)  is  the  constant  of  integration.  Substituting 


the  uniformly  convergent  series  (2.7)  into  the  integral 
yields 
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/  -^drp  =  l  J  (^)drp  Pn(cos<IO 


n=0 


(2.103) 


=  -  l  ~r  Pn(cosil>)  +  r  cos'!'in^.  +  b(r  ) 

n=0  n_1  rn_1  n  s  rs  s 


n#l 


where  r  >  rt  .  Bence 

p  5 


I  — P  (cos'!)  =  -i-  r  cos^in  2(rp  +  ^  -  ^coe'P) 
n-1  n-1  n  s  r 

p  v  (2.104) 


n=0 

n*l 


+  b(rg)  -  a(rg) 


This  holds  for  any  r?  >  rs#  in  particular  as  rr-*e»,  so  that 


lim(  £  -  r  )  +  0  =  -r  cos<M-n4  +  b(r  )  -  a(r  )  (2.105) 

oc  V  °  S 


he writing  i  -r^  as 


l  -  r. 


A .  ajj.  -  ?r|Be<>a«.  -  i 


(2.106) 


l*Hopital's  rule  gives 


lim  <f  -  r  )  »  -rg  cos  (2.107) 

rp-*. » 

Therefore  the  total  constant  of  integration  is  given  by 
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b(rs)-a(rg)  =  rgcos^an4-l)  (2.100) 


Putting  (2.108)  into  (2.104)  finally  yields 


P  (COSty)  ■ 


(2.109) 


-  Jt  -  rscos in 


V  l  -  rs  cos4> 


rg  cos  <fi 


Siailarly,  the  uniformly  convergent  series  (2.8) 
substituted  into  the  integral  of  eguation  (2.102)  results  in 


/  ^  drp  »  l  I  (^E.)n+1  drp  P  (cos  1 0 
*  n=0  s  n 

»  n+2 

"  L  575-  ^n+l-V n(cos®  +  c(rs) 
n=°  rg 


(2.110) 


where  c(rt)  is  the  integration  constant,  and  r  <  r,  .  with 
(2.102)  this  becoaes 


n+2 

J  -ly  rPVj—  pp(cos  ip)  -  l  +  rs  cos<J/in2(rp+  l  -  rgcos\0 


n=0 


(2.111) 


+  a(r  )  -  c(r  ) 
s  s 


This  also  holds  for  r,*r,  if  v  f  0,  since  then 

P„  (cos  y  )/(n+2)  *  0(n*Vl)  and  the  series  converges.  Hence 


rs  l  T72*a(cos®m  ^  +  rscos^£n2(rs+£s-rscos«p) 
n-0 


♦  a(rg)  -  c(rs) 


(2.112) 
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where  i,  =r,V 2-2cos v .  Bat  the  first  two  teras  on  the  right 
side  of  (2.112)  can  be  expressed  as  a  series  according  to 
(2.109)  with  rp  =rs  (yfO).  Thus 


a(rs)  -C(rs) 


=  rs  L  STF  V®08*0 

1  (2.113) 

+  rs  n^l  pn(cos  0  +  rs  cos^-  rs  cosi|,tn4rs 
n#l 


Patting  this  into  (2.111)  and  coabining  the  series ,  we 
obtain 

oo  rn+2 

J0  IVnTT  +  n=I>-  5?3  Pn(cos*>  = 

n*l  s 

3 

=  -  4,  -  rscos<|>  £n  +2r~  -  ( -|-rs  “  "S'  *  cos  (2  - 1 14> 

*rP  s 

r 

+  rs  cosi/»  £n  jJi- 


Hecalling  the  definition  (2.85)  of  the  kernel  E  and  using 
(2.109)  and  (2.114),  we  finally  arrive  at  the  desired  series 
re  presen tat ion ; 


l 

n=0 

n^l 


3L  £s  p  (cocilO  +  Ts  c.os^ 

n-1  n+1  ynuosv; 

rp  P 


*  rp  ”  rs 


E  =  < 


»  rn 

t  r  rs  (  3:  +  ^  I  p  ( cos  tlO  +  (2.115) 

L,1  rn  ^n+2  n^T'  T+5  n+1  J  *Vcosw  + 


n*0  ‘P 
n#l 


l- 14 - t <£>’ ♦ 1” % 1  •  vr«  • 


*  *  0 
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The  kernel  E ,  as  a  function  of  the  point  P,  is  haraonic 
everywhere  except  on  the  line  segaent  rp  &  r , ,  V  =0  where  it 
is  not  defined  and  where  in  any  case  we  will  have  no  need  to 
evaluate  it. 

Substituting  the  series  (2.115)  into  the  foraula  for  the 
disturbing  potential  (2.84), 


■  I  l  t  ?  (6  ,X  )  ,  r  >  R  (2.116) 

rn  nra  nm  p  p  p 


n=0  m=-n  P 


where 


nm 


(  rr25TI  Jj  Vn-"iX-'  ^  Ynm(e’X)do’  n  *  1 

ar*  t  //V(0’X)  if-?im<9'X)do  *  n  =  1 


(2.117) 


and  for  points  below  the  bounding  sphere: 


n 


TP  ‘  Jo  J- .b“<rP>  '-‘VV  1  rP  '  B  '2-"8’ 


where 


bnm(rp)  =  ?F2^T  *n^T  [}  vZn+T_  Ynm(0’X)da 


rP 


+  //  "iff  <ot  +  i£r>  - 

n  *  n  r 


n  t  1 


bim(rp>  ‘  ?lo<e>x)d<’ 


*  //  t4  -4  <£>*  +  *»  j£  >  7.m<9"»d'’' 

0  O  ** 


n 


(2.119) 
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with  the  spherical  approximations  as  defined  for  (2.56)  the 
expansions  of  the  height  anomaly  follow  immediately  from 
(2.116)  and  (2.118):  for  r,  >  B, 


1  ^  (:r-)n-1  z  7  (9  > 

~  nm  n m '  n  *  r* '  * 


s-  *  i  <r- 

n®0  m«-n  P 


nm  nn  p’  p'  *  znm  “  %  tnm  <2*120) 


and  for  r^<  B, 


ip  ‘  Jo  mi„!'""'rP)  Tnm<VV  •  y„n<rp)=S-bnm(r  K2.121) 


The  above  series  and  the  corresponding  series  (2.97)  for  the 
gravity  anomaly,  although  formulated  specifically  for  points 
below  the  bounded  sphere,  nay  be  regarded  as  convergent 
series  for  the  height,  or  gravity,  anomaly  anywhere  in  the 
exterior  spaae,  since  they  revert  to  (2.120)  and  (2.95)  for 
rp  >  R;  that  is  , 

y  (r  )  »  (— )n“a  z  ;  c  ( r  )  -  (  — )n+2  g  ; 
ynmv  p  vrp'  nm  *  nrav  p  vrp'  Bnm  ^ 

rp  >  « 

finally,  the  downward  continuation  error  of  the  height 
anomaly  is  obtained  by  subtracting  (2.121)  from  the  series 
(2.120)  truncated  at  n: 


e<Cn>  =  I  S  e  ( r  )  Y  <6  X  >  + 
p  n=0  m— n  nm  p  nra  p  p 


+  7  j  e'  (r  )  Y  ( 8,  X) 

n=ff+l  m=-n  nn  p  nm  p  p 


(2.123) 


where 


enm(rp)  *  (^)n_1  iSr^m  -  Tfr  bnm<rp)  = 


(2.124) 


65 


/ 


*  rpR2  l 
~  4tkM  2n+l 


and 


//  ^KT2<^>n+1+n4r^>n-^^„4r 

y„m(B,X)dci  .n^l 


rr>l 


e'nm<rp)  *  T&  bnm(rp)-enm(rp 


*  [JflJ  \  V  ,  ^  /  V*V-*  ,  **  T-  A. 

'//  +?J-  tn  1  Jim<e.»dc  . 

M  —  1 

(2.125) 


n  =  1 


'  ,  .  K"  ,  K  v 

(rp)_  <M(r  5 


R1  /  R  xn-lt 


2.3.2  The  Description  of  the  Hodel 

is  with  the  volumetric  density  aodel  certain 
simplifications  of  the  earth's  surface  are  required  in  order 
to  evaluate  the  integrations;  however,  this  model  places 
alaost  no  Units  on  the  conplexity  of  the  density  function 
V  (0,A),  given  by  (2.94).  The  necessity  of  a  highly 
complex  model,  when  aiming  for  expansions  of  the  potential 
and  gravity  anoaaly  to  degree  and  order,  say  300,  is 
prescribed  essentially  by  the  Hyquist  lav.  although  this 
lav  holds  only  for  expansions  of  functions  in  Cartesian 
space,  for  the  present  purposes  it  serves  as  a  sufficient 
guide.  The  Myguist  lav  (applied  to  a  great  circle)  states 
that  the  Fourier  spectral  components  of  a  function  whose 
values  are  Known  at  a  uniform  interval  of  180°/n  can  be 
resolved  to  a  degree  no  higher  than  ii  (Bath,  1974,  p.146). 
Hence,  for  5=300,  we  should  specify  values  of  the  model  on  a 
0?6  grid.  Instead  of  equal  area  blocks,  as  for  the 
voluaetric  density  aodel,  equiangular  blocks  delimited  by 
coordinate  lines  will  be  used  here  so  as  to  take  advantage 
of  the  Fast  Fourier  Transform  of  data  along  bands  of 
constant  latitude  (see  below)  and  thereby  to  make  the 
computations  manageable.  However,  the  convergence  of  the 
meridians  toward  the  poles  is  accompanied  by  an  increase  in 
the  concentration  of  the  data.  This  nonunif ormity  of  the 
data  implies  a  somewhat  larger  frequency  content. 

Unfortunately,  data  sets  at  a  resolution  of  o!6»67  km 
do  not  exist,  especially  for  the  gravity  field.  Kith  the 
(180,180)  solution  of  harmonic  coefficients  (derived  by  Bapp 
from  1®x1°  mean  gravity  anomalies  obtained  in  (Bapp,  1976)) 
as  a  base,  the  coefficients  from  degree  181  to  300  could  be 
generated  from  random  numbers  which  decay  according  to  a 
specified  degree  variance  model  for  the  gravity  anomaly. 

Such  a  model  was  obtained  by  Tscberning  and  Bapp  (1974)  and, 
in  fact,  it  has  not  been  significantly  vitiated  in  light  of 
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the  sore  recent  (180(1&0)  solution  (Rapp,  1979).  However, 
it  was  decided  to  use  the  aore  coupler  "tvo-coaponent"  aodel 
(ibid.),  since  its  paraaeters  were  deterained,  in  essence, 
by  this  high  degree  solution.  On  the  other  hand,  to  obtain 
a  saootn  transition  iron  the  actual  degree  variances  to  the 
aodel,  the  (180,180)  solution  was  used  only  up  to  degree 
100: 


^(coefficients  of  the  (180,  180)  solution), 
pn 


n  =  2  ,  . 
n  =  101  , 


,  100 
.  ,  300 


(2.126) 


where  /3„  is  the  (approziaate)  frequency  response  for  the 
operator  which  averages  a  function  over  1*z1*  blocks 
(Pellinen,  1966;  V#  *0?564  in  equation  (2.63)),  the  u„w  are 
uniforaly  distributed  randoa  values  in  the  interval 
[-0.5,0.53,  b„  is  their  degree  variance, 

□ 

bn=  JnuL  (2-127) 

m=-n 


and  the  gravity  anoaaly  degree  variance  c„  is  aodeled  in 
agal*-  according  to  (see  Bapp,  1979) 


vn+2 


3 « 405(n-l )  /n  nnono^ \ n+2  ,  140 . 03 ( n-1 )  , ^ 

n+1  (0.998006)  +  (0.914232) 

(2.128) 


n  >  3 


The  division  of  the  (180,180)  coefficients,  being  spectral 
coaponents  of  the  1°  aean  anoaaly,  by  /3„  transforas  then  to 
spectral  coaponents  of  the  point  anoaaly  (see  equation 
(2.64)). 

Ue  note  that  the  coefficients  (2.126)  define  only  the 
density  distribution,  as  given  by  (2.94) ,  and  not  the 
gravity  field  of  the  aodel.  The  ezpansion  (2.94)  for  the 
density,  being  the  solution  to  the  integral  eguation  (2.91) 
if  the  earth  is  a  sphere,  obviously  does  not  correspond  to 
reality.  But  since  the  earth  is,  in  fact,  nearly  spherical 
and  a  lack  of  knowledge  necessitates  a  certain  aaount  of 
siaulation  of  the  data,  a  density  aodel  based  on  (2.94), 
though  not  optiaal,  is  with  these  arguaents  also  not 


indefensible .  (Petrovskaya  (1979)  gives  a  for aula  by  which 
the  integral  egaation  could  be  solved  iteratively,  given  the 
shape  of  the  earth*s  surface;  however,  it  is  not  necessary 
to  iapleaent  this  procedure  for  the  present  purposes.)  The 
zero  and  first  degree  terms  in  the  expansion  of  the  density 
have  been  omitted  on  the  usual  assumption  associated  with 
the  spherical  approximation  that  the  average  and  first 
moments  of  the  density  with  respect  to  the  origin  of  the 
coordinate  system  are  zero.  However,  when  the  density  is 
weighted  by  the  surface  radius  as  in  (2.117),  the  first 
moments  do  not  vanish  and  the  potential  of  the  earth  model 
does  contain  a  first  degree  term.  (Also,  a  zero  degree  term 
is  present  because  the  discrete  grid  values  of  the  density 
do  not  average  exactly  to  zero.) 

The  earth's  surface  can  be  modeled  in  the  first  place  on 
the  basis  of  1* x1°  mean  elevations,  available  as  a  global 
data  set  (provided  by  DHAAC,  1979)  and  describing  the 
surface  with  respect  to  sea  level  with  a  resolution  of 
approximately  110  km  (higher  towards  the  poles).  The 
additional  finer  detail  from  degree  1B1  to  300  again  must  be 
fabricated  by  a  random  number  generator.  This  is 
facilitated  and  improves  the  verity  of  the  model  if,  instead 
of  ascertaining  an  independent  degree  variance  model  for 
elevations,  we  invoke  the  probable  correlation  between  high 
degree  potential  coefficients  and  short  wavelength 
topography  (Lambeck,  1979,  p.590) .  This  correlation  can  be 
derived  by  assuming  that  the  high  degree  components  of  the 
disturbing  potential  are  generated  by  the  masses  of  the 
topographic  features,  including  isostatic  compensations, 
condensed  onto  a  mean  earth  sphere,  in  order  to  account  for 
the  lower  density,  /tv=1.03  g/ca* ,  of  the  oceans  relative  to 
the  crust,  "eguivalent  rock  topography"  has  been  introduced 
(Baimino  et  al.,  1973)  whereby  the  oceans  have  been 
replaced,  on  the  ocean  floor,  by  an  eguivalent  rock  layer  of 
crustal  density  ^ue  and  thickness  ,  where  d  is  the 

depth  of  the  ocean  (d  *  0)  .  The  eguivalent  rock  topography, 
H, ,  is  measured  with  respect  to  the  geoid  and  is  therefore 
negative  in  ocean  areas: 


Hr(9,X) 


H(M) 


for  land  areas 


-(1-jpOd  ,  for  ocean  areas 


(2.129) 


Lambeck  (ibid.,  p.592)  gives  the  following  relationship 
between  the  (disturbing)  potential  harmonic  coefficients  t, 
and  the  corresponding  coefficients  of  the  eguivalent  rock 
topography  h„m: 
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t 


nra 


4TTRK  yc[l-  (5zL)n] 


(2.130) 


where  /tc  is  assuaed  constant  (2.67  g/ca5),  5  is  the  lean 
earth  radios  (6371  to) ,  and  o  is  the  depth  of  coapensation 
in  the  Air y-Heiskanen  isostasy  aodel  (a  value  of  0=50  ka  was 
found  to  give  better  agreeaent  between  the  (180,160) 
solution  and  the  1*x1*  elevations  than  0=30  ka  (Rapp, 
private  coaaunication) )  . 

The  coefficients  hnm  were  determined  to  degree  180  froa 
the  1*xl*  aean  elevation  data  set  (which  includes  negative 
depths  ~d  in  ocean  areas)  and  the  definition  of  Hr  (2.129) 
and  hence  constitute  its  spectxua  for  0  £  n  *  180.  The 
spectrua  froa  degree  181  to  300,  assuaed  to  be  related  to 
the  potential  spectrua  according  to  (2.130),  is  directly 
obtained  froa  the  randoa  coefficients  (2.126): 


f  Tff  Hr<e’x)  Ynm(0’X)do  »  0  -  n  -  180 

(2.131) 

,  181  <  n  <  300 


Having  thus  obtained  its  spectrua  to  degree  300,  the 
equivalent  rock  topography  is  evaluated  on  a  0t6  by  0*6 
global  grid  using  the  expansion 


2p+i _ xtm _  ZpT 

n_1  47TK  yctl-  C^-)n]  A 


Hr(8,X) 


300 


I 

n=0 


nn 


Ynm(0’X> 


(2.132) 


and  the  actual  topographic  surface  of  the  aodel  (with 
respect  to  the  geoid)  is  then 


H( 0 , A  ) 


Hr(8,X) 

0 


if  Hr( 0 , A)  >  0 
if  Hr(  0  ,  A  )  <  0 


(2.133) 


The  degree  variances  of  hnm  are  shown  in  Pig.  13.  Although 
the  coefficients  h**,  ,  Of  n*  180  refer  to  aean  topography,  no 
unsaoothing  of  this  portion  of  the  spectrua  was  perforaed. 
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as  with  the  potential  spectra#  (see  (2.126)),  since  it  would 
hawe  destroyed  the  essentially  continuous  transition  to  the 
■odeled  high  degree  part  of  the  spectrum.  These  types  of 
manipulations,  while  perhaps  not  strictly  acceptable  from 
the  theoretical  viewpoint,  are  designed  to  produce  a  model 
which  renders  as  faithfully  a  representation  of  the  earth  as 
possible  including  the  requirement  that  the  spectra  of 
gravity  and  topography  decay  with  no  major  jump 
discontinuities . 

The  coefficients  g„w  (eguation  (2.126))  find  further 
utility  in  the  definition  of  the  geoid  height  N(0,X) 


b  3  300  n  ~ 
n=2  m=-n 


This,  of  course,  is  not  the  true  undulation  of  the  earth 
model,  since  the  set  of  coefficients  is  not  the  gravity 

anomaly  spectrum  of  the  model;  but  (2.134)  serves  well  * 
enouah  for  the  present  purposes.  Pinally,  the  geocentric 
distance  to  the  model  surface  is  given  by 


r  ( 0 , X)  =  p( 9 )  +  N( 9 , X)  +  H( 9 , A )  (2.135) 

o 


where  p  is  the  geocentric  distance  to  the  reference 
ellipsoid  (equation  (2.81)). 

k  typical  profile  of  the  resulting  model  surface  is 
depicted  in  Pig.  14,  and  contrasted  with  the  5°  and  1°  mean 
values  of  r$ .  The  farthest  distance  of  the  modeled  surface 
from  the  origin  of  coordinates  is  rs  (mar) =6361989. 115  m,  at 
a  latitude  of  -9l3.  The  radius  of  the  bounding  sphere  was 
rounded  to  B=6382000  a.  The  BBS  deviations  of  the  surface 
from  this  bounding  sphere  range  from  3760  m  at  the  equator 
to  25100  m  near  the  north  pole. 


2.3.3  The  Equations  for  t&e  Numerical  Analysis 

The  integration  of  the  density  integrals  is  tractable 
only  if  r,  and  y  are  assumed  constant  within  each  0?6x0t6 
eguiangular  compartment.  The  resulting  salient  surface  does 
not  satisfy  Liapunov's  condition  of  continuous  curvature 
(Gunter,  1967)  and  the  convergence  of  the  potential  and  its 
first  derivatives  to  their  respective  values  on  the  surface 
is  not  guaranteed.  However,  the  need  for  evaluations 
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exactly  on  the  surface  will  generally  not  arise  in  these 
numerical  studies.  We  then  have,  for  example,  froa  eguation 
(2.117) 


*nm 

N-l 

'  17-sfei  Jo 

2N-1 

r 

L 

J=o 

•  w  t- 

R 

0i+1  /j+1  cos  m  X 

J0  Pn|m|(cos0)  sin9d0/x.  (sin|mi 

i  J 

x)dX 

(2. 

136) 

where 

rjK  =1/(n-1) , 

n  >,  0 , 

n  f  1,  and  -vj,  =1; 

and  where 

8i 

=  i A8  ,  i  = 

0,  1, 

.  .  .  ,  N  -  1 

N  - 

180° 

9 

—  (2- 

.137) 

Xi 

=  JdA  ,  j  » 

0,  1, 

...»  2N  -  1 

where  =0*6=  nX  ,  N=300.  Sote  that  rs  and  v  are 

evaluated  at  (0;., x  ,  )  ,  i.e.  at  the  center  of  each 

coapartaent.  *  * 

The  coaputation  of  all  (N+1)r  coefficients  tnw,  is 
virtually  impossible  for  large  N ,  such  as  N=300,  without  the 
aid  of  the  Past  Fourier  Transform  (FFT) .  This  is  an 
extrenely  efficient  routine  for  computing  the  Pourier 
spectrum  of  a  function.  In  our  case,  the  PPT  is  applied  N 
tines,  once  for  each  of  the  latitudinal  bands  9t -constant 
in  which  r3  and  v  are  functions  of  longitude  only,  host  of 
the  following  derivations  are  found  also  in  a  aore 
generalized  fora  in  (Coloabo,  1961)  ,  but  are  included  here 
for  the  sate  of  completeness.  Denote  by  X^^ithe  integrals 
of  the  Legendre  functions: 


Wl  -  Je  Pnm(cos9)  sine  d0  (2.138) 


These  can  be  computed  using  the  recursive  algorithm  of  Paul 
(1978)  and  by  noting  that 

=  ^mi  ,  i  -  0,  1,  .  .  .  ,  JN-  1  (2.139) 


It  has  oeen  found  that  Paul's  recursives  are  sufficiently 
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accurate  (for  the  tests  conducted  here)  up  to  degree  300  for 
all  latitude  intervals  (Colombo,  personnal  communication) . 
Mow  let 


F  .  .  * 
nij 

‘  vCOi 

+*’ 

xj  +i) 

f  Ts@i+i 

1  R 

*  X  i+£) 

(2. 

140) 

then 

we  have 

Rnn 

N-l 

2N-: 

L  ,\l+l _ , 

^nm 

4n(2n+l) 

.  I  x 

i=0 

n  H 

i  I 

j=0 

F^!XJ 

/  CU2>  111  A  .  .  * 

(sin!ra!X)dX 

(2. 

,14  1) 

And 

if  we  let 

Xj+1 

cos  m  XdX 

amj  = 

cos  m 

*  amj 

“S 

bmj 

sin  ra 

XJ 

«  fAJ+1 

sin  m  XdX 

(2. 

.142) 

then  using  (2.137)  and  the  angle  sun  formula  for  sinusoids, 
we  f  ind 


q  .  ”  B  el  •  —  A  j 
mj  ra  mj  ra  raj 


emj  =  Am  araj  +  Bm  braj 


(2.143) 


where 


m 


m  > 

0 

• 

f 

1 

13  =  i 

— b  1 
ra  ml 

9 

V 

o 

m  - 

0 

IQ  J 

AX 

,  m  »  0 

(2.144) 


low  let 
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P_mS  2N-1 

fq  »'*  X 

4nnu  j»0 


_  rcos  m  X  -U 

rnij  ^sinmXj 


(2.145) 


then  with  (2.142)  and  (2.143),  (2.141)  becomes 


Rr)  --  — 

4i<2n+l)  i^Q  xnmi  ^Bm  Pnmi  "  ^  qnmi^  *  m  ~ 


N-l 


nm 


Rnn 


N-l 


(2.146) 


^.-m  =  4ir^2n+l)  i|Q  xnmi  ^Am  Pnmi  +  Bm  qnmi^*  m  *  0 


The  Fooriec  spectra  (2.145)  are  coaputed  siaultaneously  by 
the  FFT;  but  only  relatively  minor  savings  in  tiae  are 
achieved  with  respect  to  the  standard  aidpoint  nuaerical 
integrations,  because,  as  eguations  (2.145)  stand,  the  FFT 
aust  be  applied  once  for  each  i  and  n,  i.e.  a  total  of  H1- 
tiaes.  A  substantial  redaction  in  the  conputational  tiae 
results  by  introducing  the  binoaial  expansion  of  the 
functions  Fwi^.  Letting 


h(e,x)  =  1 


rs( 9,1) 
R 


(2.147) 


ve  have  the  usual  binoaial  series 


n 

(1  -h)n 


□ 


(“)  (-l)k  hk 


(2.148) 


Since  0<h<<1,  it  is  not  necessary  to  take  aore  than  a  few 
teras,  say  K+1,  in  order  to  achieve  sufficient  accuracy. 

The  Pourier  transfer  as  pBW>i,  then  becoae 


{ 


nmi, 


‘nmi 


K 

k«0 


(-1)J 


„  2N-1 

<£)  l  * 

*  j-0 


kij 


r  cos  m  X  u , 
‘ sin  m  \y 


(2.149) 


where 
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Fkij  “  v(9i+i  »  Xj+iJ  [h<6i+* 


(2.150) 


requiring  a  total  of  (K+1)N  applications  of  the  FPT. 

The  density  integrals  corresponding  to  the  downward 
continuation  error  in  the  gravity  anonaly  and  height  anonaly 
are  treated  similarly.  The  factor  depending  on  the  powers 
of  the  surface  radius  in  the  case  of  gravity  anoaalies 
(eguation  (2.100))  is  expanded  as  follows,  where 


h(0,X)  =  1 


(2.151) 


Then 


<~-)n  -  <£fi)n+1  *  (1  -  h)~n  -  (1  -  h)n+1  = 
rp  rs 

°  (2n+l)  [E  +  -  (53  +  E1*)  +  S5  +  .  ..) 


we  note  that  the  series  (2.153)  represents  also  the  first 
degree  tern  in  (2.124). 

If  we  substitute  the  series  in  h  (2.152)  into  the  error 
coefficients  d*w(r?)  (eguation  (2.100)),  then 
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dnm(rp)  "  *£r/(  V<0*X)  fi  *nm(0’X)do  +  aQm(rp)  (2.155) 


where  d„m(r»)  represents  the  part  of  the  error  coefficient 
containing  n  to  the  third  and  higher  powers: 


d 


nm 


//  v(6,X)  [Sl®±li  <hJ 

02 


+  E“ )  + 


1  • 


(2.156) 


•?„^0,X)da 


The  total  downward  continuation  error  (2.99)  becoaes 


<4’8)  =  ?  jc  J-n  »■  R  5n»(e’A)da  ?™<VV 


+  „io  J-n  3-<rP)  VV  + 


+  l  I  a-  (r  )  Y(6  .AJ 
n«=n+l  m=-n  ntn  P  nm  P  p 


(2.157) 


where  d^(rr )=dw„  (rp) -(B/rr)"4t  gK„  .  The  value  of  the  first 
tern  above  is  zero  if  the  point  P  lies  on  or  above  the 
earth's  surface.  To  show  this,  let  rF  be  fixed  (thus  fixing 
the  set  art  )  and  consider  the  function 


D(9,X) 


v<0,X)E 

i 

0 


(0,1)  e  02 

(0,X)  e  Oi 


(2.158) 


(recall  that  h*  1-rf /r,  (  \  )).  The  function  D  is 

continuous  and  is  therefore  expandable  as  a  series  of 
spherical  ha  monies,  with  a  change  of  notation,  we  can 
write 


D(W  "  nL  J  W  U  V<0*X>  E  5nm(0'x)i»* 
n*u  id*— n  02 

•  Y  C0,!)  (2.159) 

nm  p  p 
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where  ( 8P, Ap)  is  any  point  on  the  sphere  of  radius  rp ,  and 
where  the  coefficients  of  this  series  are  deterained  by 
substituting  (2.15S)  into  (1.6).  But  since  the  error  series 
is  to  be  evaluated  only  for  points  in  er,  ,  the  series 
(2.159),  hence  the  first  ter*  in  (2.157),  vanishes.  When 
foraulating  the  downward  continuation  error  for  the  aean 
anoaaly,  the  saae  arguaents  apply.  In  this  case,  we 
consider  the  average  of  the  function  D  over  a  spherical  cap 
on  the  sphere  of  radius  rp  (the  set  cr,  thereby  decreases  in 
aeasure)  . 

A  siailar  reasoning  clearly  leads  to  a  downward 
continuation  error  in  the  height  anoaaly,  where  the 
quadratic  tern  hL  is  absent: 


e<S„>  -  l.  I  hm<rp>  ?„o,<W 


n=0  m=-n 


+  I  l  5 '  ( r)  f  (0  .)) 

__ - . .  _  nra  P  nm  p  p 

n=n+l  m=-n 


(2.160) 


where 


e  (r) 
nm  p 


02  P 

,?nm  (9’A)do 


(2.161) 


and 


5Dm<V  '  5„„,<V  -  TS  (I;)"'1  t„m 


(2.162) 


Eguations  (2.156)  and  (2.161)  are  essentially  equivalent  to 
those  of  Petrovskaya  (1979)  -  if  we  asuae  rt*rp»1,  then 

e  (r  ).  =  -1)  e  (r  )i 

nm  p  I  equation <2. 161)  (n+4?  nm  p  ‘Petrovskaya  (2*163) 
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Since  the  function  D,  defined  by  (2.158),  has 
discontinuous  derivatives  at  the  curve  separating  the 
surfaces  cr,  and  o*v,  the  series  (2.159)  converges  very  slovly 
to  zero  for  points  near  the  earth's  surface  (a  case  of  the 
well  known  Gibbs'  phenomenon) .  This  seens  to  be  the  reason 
for  the  increase  by  one,  two,  or  nore  orders  of  nagnitude  in 
the  sun  of  the  first  5  error  coefficients  (2.100)  (or 
(2.124))  when  the  linear  (or  quadratic)  tern  in  the  h-series 
is  retained. 

With  the  assuned  constancy  of  rt  and  v  over  the  0?6  by 
0?6  compartment s,  equations  (2.156)  turn  into 


.2  N-l 


^nn/rp^  475*1  xnmi  tBm**nmi  "  Am  ^nmi^  ’  m  -  0 


)  2  N-l 


(2.164) 


^n,-n/rp^  47rr2  xnmi  ^Am^nmi  +  Bm  ^nmi^  ’  m  >  0 


where  the  Fourier  transforas  are  given  by 


2N-1 


<r>  -  x  i=o 

Hnmx  k=3  j-0 


2" A  ,  =  rcosml-ii 

^  Aij  Fkij  ‘  sin  niA.' 


(2.165) 


and  where 


Fkij  -  ^i+i’W  lB<8i+*#^+*>3  * 


_  _  n(n+l)  „  _  n(n+l)(n2+n+18) 

xn3  xn4  S  »  xn5  - 125 - 


(2.166) 


and  finally,  where 

1  ,  if  rp  <  rs  (®i+i> 
0  ,  if  rp  i  rs  <ei+$> 


(2.167) 
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The  factor  Aij.  in  (2.16  5)  ensures  that  only  points  on  the 
surface  above  the  sphere  of  computation  (radius  rp) 
contribute  to  the  error*  and  only  those  latitudinal  bands  in 
which  at  least  one  such  point  exists  are  included  in  the 
sums  (2.164).  The  coefficients  of  the  downward  continuation 
error  in  height  anomalies*  equation  (2.161)*  then  also 
become 


R2rp 


N-l 


I  xnm,  [Bn  pnmi  -  Am  a  ,  ]  ,  m  2  0 


'n,m  8itkM  Anmi  1  n  *nrai  m  ^nmi 


R*r 


N-l 


Jn,-m  8irK®  xnmi  Pnmi  +  Bm  ^nmi^*  m  >  0  (2.168) 


where 


I  pnmi  i 
^nmi 


K  2N-1 

I.  snk  *  4 

k=4  j=0 


rCOSmXji 

ij  *kij  ‘sin  m  X  j  r 


(2.169) 


and  where 


Fkij  =  FJ  rs<9i+i’  Aj+i5  Fkij 


n(n+l)  ..  _  n(n+l)  _  n(n+l)(n2+n+18) 

yn4  =  ~T5  *  yn5 - 15“’  yn6  - 555 - '  ** 


(2.170) 


2.3.4  Xhe  Numerical  Besolts  of  fhe  Error  Analysis 

JLs  in  the  case  of  the  volumetric  density  model*  the 
radial  dependence  of  the  error  "coefficients'*  (2.100)  and 
(2.124)  precludes  the  mapping  of  the  error  on  a  global 
scale,  on  the  other  hand*  because  the  earth's  surface  is 
more  or  less  rotationally  symmetric*  a  reasonably  thorough 
analysis  lies  within  feasibility  if  we  limit  the  latitudinal 
extent  of  our  investigation  to  a  few  representative  degrees. 
Table  1  lists  the  regions  to  be  considered*  as  well  as 
several  representative  radii  tor  each  range  of  latitudes* 
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Table  1:  The  Region  of  the  surface  aodel  in  which  the 
downward  continuation  errors  are  computed. 


No.  of  points  on  all 
spheres  of  radius  R, 
such  that  0^rp-rs<  t 

\ 

300  b 

t  -  100  B  » 

50  B 

l  10  *  j 


North  Lat.  South  Lat. 


Region* 
defined  by 
latitudinal 
range  of 

Represen¬ 
tative 
radii,  H, 

w 

Longitudinal 
ranges  of 

1 

6374550 
6374  850 
6375150 

[0“.3,3594.7] 

12  th  i  |y;l  £  24?3 

6  375450 

6  375750 

6  376050 
6376350 

[69°. 9,  89t7] 

6  376650 
6376950 
6377250 

^90tl,309t9] 

11 

6  363350 

6  363750 
6364150 

[o!3,359°.7] 

47“  1  i  |f(.| *57^3 

6364550 
6364950 
6365  3  50 

6  365750 

[60°.  3,  80  !l] 

6  366150 

6  366  550 

[2  60 1 1,28  5  “3] 

111 

6  357250 

6  357550 
6357850 

[o#.3,359*.7] 

75^3  £1^1684! 3 

6358150 

6  358  4  50 

[l70°.1,260tl] 

683  714 

234  273 

87  161 

17  0 


North  Lat.  South  Lat. 


North  Lat. I  South  Lat. 


2090 

572 

171 

19 


Pigure  15:  Subregions  over  which  RBS  downward  continuation 
errois  are  computed. 
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see  also  Fig.  15.  He  note  that  foe  each  constant  rf  =Br  ,  the 
coefficients  (2.100),  or  (2.124),  correspond  to  the  error  at 
any  point  which  lies  on  the  sphere  of  radius  B r  and  above 
the  earth's  surface.  Therefore,  instead  of  evaluating  the 
error  series  at  a  few  isolated  points  on  the  surface,  a  sore 
characteristic,  although  less  realistic,  assessment  is 
obtained  if  the  series  is  evaluated  on  the  given  spheres  at 
all  those  points  which  are  removed  from  the  surface  by  no 
more  than,  say,  100  ■  (see  Pig.  17).  Table  1  also  shows  the 
numoer  of  points  on  a  G?6x0?6  grid  that  lie  on  the  spheres 
having  radii  and  within  a  tolerance  of  300  a,  100  a,  50 
a,  or  10  a  above  the  modeled  surface. 

The  aaount  of  calculations  can  be  further  reduced  by 
requiring  no  aore  than  two-  or  three-digit  accuracy  in  the 
error  estimates;  any  atteapt  at  greater  coaputational 
accuracy  is  unavailing  and  therefore  unjustified.  Hence,  we 
aay_accordingly  restrict  the  number  of  teras  of  the  series 
in  h,  (2.152)  and  (2.153).  Table  2  correlates  the  accuracy 
of  the  series  with  the  nuaber  of  included  teras  for  the 
worst  cases  in  each  of  the  regions  of  Table  1.  The  worst 
case  oocurs  when  n=500  and  h=max  (h{.)  .  Similarly,  the  series 
in  h  (2.148)  may  be  limited  to  several  teras  (see  Table  2). 
Numerical  tests  confirmed  that  the  relative  accuracy 
guaranteed  for  the  truncated  series  in  h  was  not  degraded  in 
the  process  of  determining  the  error  coefficients. 

Equation  (2.146),  with  the  approxiaate  transforms 
(2.149)  substituted,  was  used  to  determine  the  spectrum  of 
the  gravity  field  model  on  the  bounding  sphere  (radius 
Rs6382000  a) .  The  series  in  n  was  developed  to  K=10  thus 
ensuring  six -digit  accuracy.  The  corresponding  degree 
variances  are  shown  in  Pig.  16.  The  series  (2.95)  and 
(2.120)  for  the  gravity  and  height  anoaalies  above  the 
bounding  sphere  then  provide,  in  conjunction  with  the  error 
series,  the  corresponding  inner  series,  which  converge  below 
the  bounding  sphere,  naaely  (2.97)  and  (2.121).  All  series 
were  truncated  at  degree  ?00. 

Once  again,  it  aust  be  remembered  that  individual  teras 
of  the  error  series  do  not  indicate  errors  in  the  harmonic 
constituents  of  the  potential  (or  gravity  anoaaly),  since 
the  inner  jeries  is  not  a  spectral  representation  of  the 
field  on  the  surface.  In  the  strictest  sense,  the 
difference  between  truncated  inner  and  outer  series  is  just 
as  aeaningless,  since  it  implies  a  comparison  of  "bands*  of 
frequencies.  The  comparison  is  valid  only  if  the  inner 
series  has  converged  with  sufficient  accuracy  to  the  true 
value  being  estimated.  Otherwise,  when  limiting  the 
evaluation  of  the  downward  continuation  error  series  to 
terms  of  degree  no  greater  than  5,  we  aust  be  aware  of  the 
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Table  2:  Accuracy  of.  series  (2.152)  (without  the  linear 
ten)  (  series  (2.153)  (without  the  quadratic 
tern) ,  and  series  (2. 148)  versus  the  nuaber 
of  included  terns. 
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effect  of  truncation,  which  is,  of  coarse,  not  known,  with 
the  present  density  model,  it  appears  likely  that  the  usual 
degree  variance  aodels  of  the  earth *s  gravity  field  will 
provide  estimates  of  the  neglect  of  higher  degree  terms, 
although  this  is  not  evident  frow  equation  (2.101)  or 
(2.125).  If  n  =  300 ,  then  the  BUS  (root  wean  square)  values 
of  the  truncation  error,  based  on  the  aodel  (2.128)  and  a 
spherical  earth,  are  about  36  ca  for  the  height  anoaaly  and 
a  considerably  sore  significant  30  agal  for  the  gravity 
anoaaly.  The  root  aean  squares  of  the  evaluations  of  the 
error  series  truncated  at  n=300  over  each  of  the  regions  of 
Table  1  are  shown  in  Table  3  for  the  gravity  anoaaly  and  in 
Table  4  for  the  height  anoaaly.  Also  shown  for  coaparison 
are  the  BflS  values  of  the  truncated  series  of  the  respective 
anoaalies  theaselves.  Some  aaxiaua  absolute  values  of  the 
error  are  also  shown  parenthetically.  Note  that  the  units 
in  the  error  columns  are  /tgal=10“*  m/sl  (Table  3)  and 
yW.a=10'6a  (Table  4). 

A  perusal  of  these  tables  indicates  that  the  errors  are 
generally  insensitive  to  the  distance  t  of  the  evaluation 
point  from  the  aodel  surface,  if  it  is  300  a  or  less.  The 
BhS  values  of  the  anoaalies,  on  the  other  hand,  show  a 
(disconcertingly  strong  and  unexplained)  decrease  as  the 
point  of  evaluation  approaches  the  surface.  There  is,  of 
course,  the  expected  increase  in  error  with  increasing 
latitude;  and  within  each  latitudinal  range,  the  BBS  errors 
show  a  significant  increase  (not  shown)  with  each  decrease 
in  radial  distance  B P .  Soae  correlation  between  the  error 
and  large-scale  topographic  features  can  be  detected.  Por 
example,  the  error  over  central  South  America  is  generally 
twice  as  large  as  in  the  topographically  lower  Caribbean 
Sea.  Here,  in  turn,  the  errors  are  almost  double  those  over 
both  the  Indian  Ocean  and  subcontinent  (just  south  of  the 
Hiailayan  massif),  showing  that  the  correlation  is  somewhat 
inscrutable.  Wo  definite  correlations  with  topography  are 
discernible  (by  visual  inspection)  from  the  values  listed 
for  the  subregions  in  the  midlatitudes  and  polar  areas. 

Table  5  shows  how  the  errors  accumulate  in  steps  of  30 
degrees.  The  accumulation  is  generally  monotonic  for 
n£  270,  but  a  noticeable  downward  trend  occurs  at  n*30Q.  In 
view  of  Pigures  6  and  7  this  is  not  an  unexpected  feature. 

To  eliminate  the  difficulty  of  interpreting  these 
results,  we  could  consider  a  smoothed  gravity  field,  where 
the  resolution,  by  definition,  is  limited  to  the  first  300 
harmonic  constituents.  In  order  to  make  the  meaning  of 
limited  resolution  tor  the  inner  and  outer  series 
comparable,  we  identify  the  cutting  oft,  or  filtering,  of 
higher  frequencies  with  a  weighted  averaging  process  as 
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Table  3:  BBS  downward  continuation  error  of  spherical 

harsonic  series  of  gravity  anomaly  and  BBS  gravity 
anoaaly  evaluated  (first  300  terss)  at  the  points 
of  Table  1  and  at  the  indicated  subregions. 

Baxiaua  absolute  values  for  each  region  and  sub- 
region  are  given  parenthetically. 


Region 


6  9  °  9- 

89°.7 


.28  .26 

.33  (2.0)  .26  (1.0) 
.35  .29 

.33  (2-0)  .33  (.73) 


.46  .27 

.48  (6.1)  .33  (1.1) 
-51  .31 

.47  (4.1)  - 


60 1 3- 


359.7 


80  "  1 


21. 

22.  (57.) 
22. 

29.  (57.) 


BBS  ag 

0  £  n  £  300 
[■gal) 


290t 

1- 

309°.  9 

.41 

.50 

.43 

(1.8) 

.30 

(-30) 

1.1 

1.1 

1.2 

(4.4) 

1.1 

(4.1) 

260*. 

285*.3 

16. 

- 

16. 

1b. 

(41.) 

16. 

(33.) 

19. 

21. 

21. 

(46.) 

23. 

(46.) 

25.23 


24.36 


290 ! 1- 
309*.9 

43 

.38 

44 

.38 

32 

.14 

102 

.98 

27 

.27 

28 

.59 

26 

.73 

28 

.18 

60.3- 
359  *.7  !  80  !l 


28.4 

28.6 

29.2 

27.8 


23.4 
23.3 

23.5 
26.2 


20  .4 
18.3 
18.9 
9.7 


31.3 

33.9 

34.6 

37.0 


33.4 

31.6 

35.6 
37.2 


19.4 
19.7 

21.5 
19.4 


am 

170ll  -  260l  1 

0^3-359* 7 

86. 

85. 

22.08 

25.17 

82.  (290) 

81.  (240) 

22.24 

24.36 

81. 

77. 

22.16 

25.67 

82.  (290) 

82.  (240) 

20.78 

23.60 

66. 

88. 

25.03 

28.39 

86.  (260) 

86.  (260) 

24.10 

28.86 

84. 

87. 

22.48 

27.81 

85.  (210) 

88.  (190) 

21.05 

16.85 

Table  4:  BUS  downward  continuation  error  of  spherical 

haraonic  series  of  height  anoaaly  and  R!1S  height 
anomaly  evaluated  (first  300  terns)  at  the  points 
of  Table  1  and  at  the  indicated  sabregions. 
daxinua  absolute  values  for  each  region  and  sub- 
region  are  given  parenthetically. 


eg  ion 


ihm 


(■ax  £(  5  ) ) 
G<  n  6  300 


0  <  n  &  300 


0*-3-  « 

359.7 

69^9* 

89*7 

290tl- 

309t9 

.26 

.23 

.43 

.32  (2.1) 

.24 

(.93) 

.53  (1.8) 

.33 

.28 

.45 

.30  (1.5) 

.36 

(-88) 

.28  (.28) 

.51 

.24 

1.1 

.53  (7.7) 

.30 

(1.3) 

1.2  (5.0) 

.58 

.27 

1.3 

.52  (5.0) 

— 

1.4  (4.6) 

o!  3- 

60°  3- 

_ 

260ll-  . 

359.7 

Bull 

285.3 

70. 

67. 

57. 

69.  (250) 

69. 

(180) 

49.  (140) 

66. 

68. 

49. 

64.  (190) 

89. 

(170) 

52.  (120) 

73. 

57. 

62. 

68.  (360) 

52. 

(150) 

64.  (140) 

63. 

47. 

62. 

71.  (360) 

46. 

(95.) 

06.  (140) 

I8BH 

’.7 

'  " 

170! 

1  -  260ll 

410 

410 

i  400  (1400) 

390  (1300) 

390 

370 

390  (1400) 

400  (1300) 

420 

420 

410  (1300) 

410  (1300) 

400 

420 

!  400  (1000) 

430  (1000) 

69 19- 

290°.  1- 

89  *7 

309*.  9 

68.12 

44.36 

67.44 

45.01 

65.73 

42.35 

62.39 

57.12 

44.89 

18.23 

40.77 

18.25 

40.12 

16.18 

— 

11.63 

6Q*.3- 

260 ! 1- 

so:i 

28513 

27.27 

37.81 

26.34 

36.81 

27.23 

37.65 

25.83 

39.13 

33.52 

33.64 

31.73 

25.19 


28.98 

27.38 

26.57 

20.96 


30.27 

32.17 

30.27 

26.94 


23.41  34.87 

24.43  35.94 

24.27  35.79 

29.32  36.86 


0l3-359*.7ll70.1-260 


14.82 

14.06 

10.64 

8.17 


30.43 

29.92 

25.39 

25.17 


5.75 

4.38 

4.20 

3.65 


41.85 

43.08 

37.62 

39.30 
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Table  5:  BBS  downward  continuation  error  of  gravity  anomaly 
and  height  anomaly  spherical  harmonic  series  in 
steps  of  30  degrees:  n«30t60 , . . . ,300,  based  on 
those  points  of  Table  1  for  which  t«10G  a. 


Region 

n 

bhs  £(ag)  >ugal 
>:  0l3  -  3591 7 

BBS  £  (  5  )  M* 

>:  0l  3  -  359*. 7 

M 

S 

N 

S 

30 

.0028 

.0079 

.0036 

.011 

60 

.015 

.034 

.018 

.046 

90 

.041 

.092 

.046 

.12 

120 

.061 

.14 

.068 

.17 

150 

.091 

.15 

.090 

.18 

1 

180 

.14 

.19 

.14 

.22 

210 

.22 

.28 

.21 

.32 

240 

.26 

.40 

.26 

.43 

270 

.40 

.50 

.39 

.55 

.33 

.48 

.32 

.53 

30 

.077 

.066 

.26 

.20 

60 

.57 

.56 

1.8 

1.9 

90 

2.1 

1.6 

6.6 

5.5 

3.4 

3.0 

11. 

9.5 

150 

5.6 

6.0 

18. 

19. 

11 

180 

9.2 

8.5 

29. 

26. 

210 

15. 

13. 

47. 

41. 

240 

Id. 

16. 

57. 

55. 

270 

23. 

24. 

74. 

75. 

300 

22. 

22. 

69. 

68. 

30 

.22 

•  ?3 

1.0 

■n 

60 

1.9 

\.r 

8.9 

90 

6.3 

4.’ 

30. 

y » 

12. 

12. 

57. 

57. 

t  ' 

21. 

20. 

100 

93. 

III 

T*  iBS 

31. 

32. 

140 

150 

46. 

46. 

220 

220 

59. 

69. 

280 

320 

270 

81. 

91. 

390 

420 

E3S 

62. 

86. 

400 

410 

69 


aatheaatically  detailed  in  section  2.2.2.  However,  the 
filter  then  cannot  perfectly  eliminate  the  effect  of 
truncation.  Because  the  Gaussian  filter  is  relatively 
efficient  in  this  respect,  it  is  chosen  here,  with 
parameters  a=15Q21,  y0  =  1S4.  The  average  is  thus  taken  over 
a  spherical  cap  of  radius  v,  *  with  values  at  the  edge  of  the 
cap  weighted  by  0.01.  The  freguency  response  (2.74)  adaits 
only  5*  of  the  haraonic  coefficients  of  degree  300  and 
thereafter  decreases  rapidly  to  zero.  The  difference 
between  the  smoothed  downward  continued  series  and  the 
smoothed  inner  series  is  then  the  downward  continuation 
error  with  relatively  little  truncation  effect.  The  BBS 
truncation  error,  based  on  aodel  (2.128)  and  a  spherical 
earth,  is  .3  agal  for  the  saoothed  gravity  anonaly  and  .7  ca 
for  the  saoothed  height  anoaaly.  The  downward  continuation 
errors  are  shown  in  Tables  6  and  7  and  follow  the  saae  basic 
pattern  as  in  Tables  3  and  4,  but  are  about  one  order  of 
magnitude  saaller  (an  obvious  conseguence  of  smoothing, 
since  the  higher  degree  haraonics,  which  are  aost  affected 
by  the  series  divergence,  are  decreased  in  magnitude)  .  Also 
included  are  the  BBS  values  of  the  saoothed  anomalies. 

Finally,  we  note  that  the  relative  downward  continuation 
error  of  the  gravity  anoaaly  is  approximately  three  orders 
of  magnitude  greater  than  the  relative  error  of  the  height 
anoaaly. 
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Table  6:  BUS  downward  continuation  error  of  spherical 

hareonic  series  of  lean  gravity  anosaly  and  BBS 
■ean  gravity  anosaly  evaluated  (first  300  terss) 
at  the  points  of  Table  1  and  at  the  indicated 
subregions.  Baxisus  absolute  values  for  each 
region  and  subregion  are  given  parenthetically. 
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Table  7:  HMS  downward  continuation  error  of  spherical 

haraonic  series  of  lean  height  anoaaly  and  RHS 
aean  height  anoaaly  evaluated  (first  300  teras) 
at  the  points  of  Table  1  and  at  the  indicated 
subregxoas.  flaxiaua  absolute  values  for  each 
region  and  subregion  are  given  parenthetically. 


3.  The  Ellipsoidal  Harmonic  Series 


There  exist  two  systeas  of  so-called  ellipsoidal 
coordinates  in  geodesy;  they  differ  in  the  definition  of  the 
latitude.  The  systea  with  the  geodetic  latitude  (defining 
the  direction  of  the  noraal  to  the  reference  ellipsoid)  is 
aost  coaaonly  used.  However,  its  three  diaensional 
generalization,  obtained  by  including  the  height  above  the 
ellipsoid  as  coordinate,  while  foraing  a  triply  orthogonal 
systea  (Bolodenskii  et  al.,  1962,  p.9) ,  does  not  yield  a 
fora  of  Laplace's  equation  that  is  solvable  by  separation  of 
variables.  In  order  to  have  any  hope  of  solving  Laplace's 
equation  witn  this  standard  aethod,  the  coordinate  systea 
aust  be  orthogonal  -  choices  of  the  third  coordinate,  such 
as  the  seaiainor  axis,  render  the  systea  nonorthogonal. 

The  spheroidal  coordinates  u,  £  ,  X  fora  a  triply 
orthogonal  systea  (Hobson,  196S,  p.421)  in  which  the  second 
coordinate,  in  geodetic  terainology,  is  the  coapleaent  of 
the  reduced  latitude;  we  aay  call  it  the  reduced  colatitude 
(in  fact,  S  i  6  ).  Although,  strictly  stated,  u,f,l  are 
oblate  spheroidal  coordinates,  we  aay  use  the-  less  precise 
noaenclature  "ellipsoidal  coordinates"  as  no  other  systea 
will  coae  under  consideration.  Their  definition  in  terns  of 
Cartesian  coordinates  is  (see  Figure  18) 


x  *  /uT+  Et  sin6  cosX 

y  »  /u 2  +  E 2  sinfi  sinX  (3.1) 

z  =  u cos  5 


where  E  is  a  paiaaeter  of  the  systea.  Froa  the  consequent 
relationship 


x2+  y  2  ,  z2  _  , 

— r  S  +  t  1 
u*fE  u 


(3.2) 


it  noted  that  the  coordinate  surface  u*constant  is  an 
ellipsoid  of  revolution  with  linear  eccentricity  E  (the 
distance  froa  the  origin  to  either  of  the  focal  points) ,  and 


92 


94 


seaiainor  axis  u.  Hence  the  set  of  coordinate  surfaces 
u^constant  is  a  set  of  confocal  ellipsoids.  The  squares  of 
the  first  and  second  eccentricities  are  given  by 


_  2  =  Ef_ 

"?+  U2 


e' 


u2 


(3.3) 


shoving  that  as  the  ellipsoid  approaches  a  sphere. 

Substituting  (3.3)  into  (2.61)  and  noting  that  al=ul+El  (see 
Fig.  18)  ,  we  find 


u  /u2  +  E  2 
/u2*  E%oss9 


(3.4) 


or  conversely, 

u2  -  i(r2-  E2)  +  E1*-  2r2E2  (1  -  2  cos28)  ]* 

Since  the  ratio  Vx* +yv /z  is  the  tangent  of  the 
spherical  colatitude  0,  equations  (3.1)  yield 


tan6  =  — ■—  -  -  tan9  (3.6) 

/u2+  E2 


Equations  (3.5)  and  (3.6)  provide  the  transformation  fros 
spherical  to  ellipsoidal  coordinates.  The  reverse 
transforaation  is,  froa  (3.1)  and  (3.6), 


r  =  /u2+  E2sin% 

tane  .  J&Z2.  tans  <3-7> 


The  choice  of  the  linear  eccentricity  was  based  on  its 
relation  to  the  seaiaajor  axis  "a"  and  flattening  f  of  any 
of  the  confocal  ellipsoids: 
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E  =  a/2f  -  f 2 


(3.8) 


For  exaaple,  ve  nay  choose 


a  =  6378140  m 
f  =  1/298.257 


giving  E -52 185 4  .4492  a. 

A  detailed  solution  of  Laplace's  egnation  in  ellipsoidal 
coordinates  aay  be  found  in  (Heiskanen  and  aoritz,  1967, 
pp41^43) .  The  solutions  for  regions  containing  the  origin 
are  P^M(iu/E)fWm (  S  ,  >  )  ,  and  for  regions  containing 
infinity,  they  are  S^Uu/B)  !„*.  (  f  ,  \  )  .  Being  interested 
priaarily  in  the  latter  region,  ve  find  the  general  solution 
to  Laplace's  equation  in  the  fora 


F(u ,  <5 ,  A  ) 


00 


I 

n=0 


♦_m  Q^|(i 
nm  npq 


E  5  ?nm(6,A) 


(3.10) 


where  the  constant  ellipsoidal  haraonic  coefficients  <pn>M 
are  necessarily  coaplez  nu abets  if  F  is  to  be  a  real 
function.  They  are  deterained  uniquely  if  the  boundary 
values  of  F  are  known  onjan  ellipsoid,  say  u=b .  Multiplying 
both  sides  of  (3.10)  by  Tnw(  S  ,  X)  and  integrating  over  the 
doaain  o’  - {(S,X)/  Oils  2^r,  0  4  S  4  n],  we  obtain  by  the 
orthogonality  of 


♦dm  '  5  !b,  Tf  II  F<f,5,X)  ?nm(8, J)  do  (3.11) 

W1  T)  0 

where  du’ssinA  d£  d>  .  m  this  way,  ve  can  define  the 
ellipsoidal  spectrua  of  F  with  respect  to  the  Legendre 
transfers.  However,  the  integral  in  (3.11)  is  not  a  surface 
integral  over  the  ellipsoid  u=b,  instead  it  is  an  integral 
over  the  unit  sphere  onto  which  the  points  of  the  ellipsoid 
with  coordinates  (  )  are  napped  according  to  the  one-to- 

one  correspondence 
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0  =  6 
A  =  A 


(3.12) 


3. 1  The  Transf oraation  f roa  Spherical  to  Ellipsoidal 

Harmonic  Coefficients 

&  brief  account  of  the  transformation  between 
ellipsoidal  and  spherical  harmonic  coefficients  has  been 
given  by  Hotine  (1969,  pp194-5) .  However,  the  formalas  are 
not  directly  amenable  to  practical  computations  beause  of  a 
lack  of  stabilizing  normalizations  in  the  formulas  and  the 
inconvenience  (from  the  geodetic  point  of  view)  of  working 
with  complex  coefficients.  A  complete  derivation  of  the 
transformation  will  be  presented  below,  which  by  introducing 
a  different  "normalization11  of  becomes  feasible,  as  well 
as  accurate,  lor  high  degree  expansions. 

The  relationship  between  spherical  and  ellipsoidal  solid 
harmonic  functions  is  established  by  using  a  more  general 
form  of  the  addition  theorem  for  Legendre  polynomials 
(Hobson,  1965,  p.364) 

P  ( vv '  -  /vz -1  /v  t 2-l  cos  w)  * 
n 


1 

2n+l 


L(-1)m  ?mn(v)  Pnm(v,)  cos  01  “ 
m=0 


(3.13) 


where  v, v'  £  C  -  (-1 ,  l]  .  If  one  or  both  of  the  variables  v,v* 
belongs  to  the  real  interval  [-1,1],  then  it  must  be 
approached  in  the  limit  through  the  complex  plane,  as  in 
(1.11).  The  formulas  (1.13)  for  Legendre  functions  with 
real  arguments  are  then  applicable.  With  v=iu/E, 
v(=cosi*iO,  equation  (3.13)  becomes  with  (1.11) 

P  (  i  cos  6  +  /l  +  3l£"  sin6  cos  u )  * 
n  t  E2 


1 

2n+l 


n 

xsio 


P__(  i  -vr)  P__(cos  5)  cos  m  u 
nm  h  nm 


(3.1«) 


Substituting  u>  =  X-t,  this  becomes  with  (3.1) 
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Pn(-^-(iZ  +  x  cos  t  +  ysint))  = 


(3-15) 


1  n  m  u 

25+T  ^ 1,11  Pnm(i  T)pnm(cos  $  cosra  <*-*> 
m=u 

If  we  aultiply  both  sides  by  cosa*t  or  sina*t  and  integrate 
with  respect  to  t  over  the  interval  [-7T,7t],  then,  since 


II 

/  cos  m ( X-t)  ( 
J  — ir 


cos  m' t 


cos  m*  X . 


dt  -  2ir  e  . 


(3.16) 


ve  obtain  eith  q= (iz+xcost+ysint)/E 

J-  f11  P  (q)  (c9smJ)  dt  =  P  (i-^-)P  (cos  6)* 

2ir  J -ir  n  H  sm  m  t  2n+l  nm  E  nm 


. cos  ra  X . 
^sin  mr 


(3.17) 


next,  consider  the  relationship  (Hobson,  1965,  p.98) 


— i-  f11  (icos9+  sine  cost)”  cos  rat  dt  *» 
ir  Jo 

n!  in-m  - -  P  ( cos  9  ) 


[(2n+l)(n+m)! (n-m)!]*  nm 


(3.18) 


We  note  that  the  integrand  above  is  an  even  function  and 
periodic  with  period  2TT  .  Bence 


TI+X 


1  n 

-7T-  f  ( i  cos  8  +  sin9  cost)  cos  m  t  dt 
2ll  TT+X 

_  ,  j  n-m  r- — 

- — _ — L£m - -  p  (cos 9) 

[ ( 2n+l) ( n+m) ! (n-m) ! ]* 


(3.19) 


for  arbitrary  X  .  How 


V- 
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^  T^X 
"2rr  “Lti+X 


(icos0+sin6  cos  t)n  cos  m  t(®?®  )  dt  = 

sin  mA 


(3.20) 


*  -k  C» <icosst  s1"8  «st)nClS!:«>d' 

+  -^-/1T+A.  (i  cos8+  sin6  cost)nsinmt  dt  (sinm^) 

^71  — TT  +  A  COS  mA 

The  second  integral  on  the  right  side  is  again  periodic,  but 
an  odd  function;  the  integral  is  therefore  zero.  Consider 
the  saae  point  (x,y,z)  as  above,  but  now  in  spherical 
coordinates: 


x  =  r  sin  0  cos  X 
y  =  r  sin  0  sin  X 
z  =  r cos  0 


(3.21) 


Then  icosO  +sin9  cos  (  X -t)  =Eg/r.  Now  uultiply  (3.19)  by 
cosa>  or  sinaX  and  substitute  (3.20) ,  thus  arriving  at 


1  f* 
"2tt  U 


n,cos  m  t.  .. 
q  (sinmt)dt 


in-m  n! 


[ (2n+l)(n-m) ! (n+m) ! ]* 
P„m(cos9)(sin^ 


(3.22) 


Mote  that  by  the  orthogonality  of  the  sinusoidal  functions, 
the  integral  in  (3.22)  is  zero  for  a  >  n.  P„  (q)  is  siaply  a 
polynoaial  in  g : 


p„(q) 


1  f  Jzlfj*.-  2jO<  v-  ,i] 

2  k=0  k!  (n-k) J (n-2k) !  4 


(3.23) 


while  q"  can  likewise  be  expressed  as  a  finite  sua  of 
Legendre  polynomials  (Hobson,  1965,  p.44) : 


v 

I 

k-0 


2n_2k(n-k)!(2n-4k+l) 
k!(2n  -  2k  +  1)! 


n-2k 


Cq) 


* 


(3.24) 


where  [x]  denotes  the  largest  integer  less  than  or  eaual  to 
x. 
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Substituting  (3.23)  into  (3.17)  and  using  (3.22),  we 
find  the  tea  ns  fori  at  ion  from  spherical  harmonics  to 
ellipsoidal  ha  monies: 


(±~) 


Ynm<6’A> 


( -1 )min  f  (2n-2k)> 
2n  k=0  <n~^) ! 


( JL)n“2k  _ ?n-2k.Tr,(9,A) _ 

E  [ (2n-4k+l)(n-2k-m) ! (n-2k+m) ! ]* 

-n  <  m  <  n  ,  s=  [— ] 


(3.25) 


And  substituting  (3.24)  into  (3.22),  sating  use  of  (3.17), 
we  obtain  the  reverse  transformation,  from  ellipsoidal  to 
spherical  hamonics: 


(TT)n  Ynm(9’X)  "  t<2n+l>  Em<n+m> ! (n-m) ! ]* 

?  2p~2^(n-k) !  =  ..  .. 

k£0  ( 2n-2k+l)  !  k!  n-2k ,  |  m|  v  1  T}  *n-2kfm*°  » A  }  (3 


26) 


-n  <  m  <  n 


s=  [■ 


n-  m 


] 


It  is  important  to  realize  that  equations  (3.25)  and  (3.26) 
represent  transformations  between  (inner)  solid  harmonic 
functions,  i.e.  they  hold  for  points  in  three-dimensional 
space.  Therefore  it  would  be  fallacious  to  deduce  the 
relationship  between  spherical  and  ellipsoidal  harmonic 
coefficients  simply  by  inserting  (3.25)  or  (3.26)  into  the 
respective  harmonic  series,  tie  recall  that,  according  to 
one  interpretation,  the  harmonic  coefficients  constitute  the 
spectrum  of  a  function  restricted  to  a  coordinate  surface, 
either  a  sphere  or  an  ellipsoid.  But  since  neither  r  nor  u 
is  constant,  respectively,  on  the  ellipsoid  or  sphere, 
(3.25),  (3.26)  do  not  provide  the  relationship  between 
surface  harmonic  functions  on  corresponding  coordinate 
surfaces. 
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To  derive  the  transf or nation  between  the  haraonic 
coefficients  of  the  two  series,  we  resort  to  their 
alternative  interpretation,  saaely  as  density  integrals. 
The  realization  of  ellipsoidal  haraonic  coefficients  as 
density  integrals  is  iaaediate  once  a  series  expansion  is 
found  for  the  reciprocal  distance.  We  have  for  uf  >  u 


1_ 

£ 


T  K  1 

n=0  m=-n 


<3.27) 


nro  p  ’  p  nm  * 


Por  u,<u,  the  roles  of  a,  and  a  are  obviously  interchanged. 
Equation  (3.27)  was  derived  by  Neumann  (1848)  and  a  much 
nore  thorough  derivation,  though  unfortunately  replete  with 
typographical  errors,  can  be  found  in  (Hobson,  1965, 
pp424— 430)  .  The  series  (3.27)  converges  uniforaly  for  up>  u 
and  can  therefore  be  integrated  tern  by  tera  when  inserted 
into  the  integral  for  the  potential: 


*  V(uD,  fiD>  Xd) 


dft  = 


fl 


X  _  ®n|m[  T?  *  Unm  *nm(6P,XP)  <3-2®> 

n=u  m=-n  1  1 


where 


u 

nm 


E(2n+1) 


a 


///  *<«,«,*>  an 


(3.29) 


Changing  to  spherical  coordinates  under  the  integral  sign 
and  using  (3.25)  yields 


,n+l 


u 


nm 


2nE 


k-0 


(2n-2k)  t  ,  R  .n-2k 
k!  (n-k) ! 

-i 


•  [(2n-4k+l)(n-2k-m)  !  (n-2k+m)  !  ]"’ k  III  P(-2r) 

«  R 


r  .n-2k 


*^n-2k,m(  9,X) 


(3.30) 
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Recalling  the  density  integral  interpretation  for  the 
spherical  haxaonic  coefficients  (2.31)  ,  this  results  in 


□m 


in+1fe^  f  ( 2n-2k) !  _  /  2n-4k+l  ~ 

2 n  k! (n-k) !  V  (n-2k-m) ! dn-2k+m) ! 


.  R  .n-2k+l 

<T> 


vn-2k,m 


(3.31) 


Both  the  Legendre  functions  and  the  coefficients  of 
the  above  sun  are  difficult  to  calculate  for  large  n  with 
computers  using  finite  digit  arithaetic.  With  the  goal  of 
aore  tractable  coaputations  in  aind,  consider  the  definition 


„n+l 

_ 


^n+1 
2nn ! 


(2n+l) 


Vi 


ia 


(2n+l)(n-m) ! (n+m) ! 


Qnm(i0P> 


(3.32) 


where  <rf  =uT/E,  <r,  =  H/E.  Then  the  ellipsoidal  series 
expansion  of  ▼  (equation  (3.28))  changes  to 


VP  '  Jo  JL  S»I»I<’p>  3n»  ?nm<W  <3-33’ 


where  now,  by  coabining  (3.31)  and  (3.32),  we  have 


unm  "  k£0  Xnmk  vn-2k,m  *  s*  (3.34) 


with 


_  (2n-2k)  !n!  f  (2n-4k+l)(n-ni)  !  (n+m)  !  1 

nmk  (2n) !  k!  (n-k) !  *■  (2n+li(n-Sk-m) !  (n-2k+m) !  2k 

o# 


n  >  0  ,  -n  s  m  s  n  ,  0  <  k  s  s 
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The  calculation  of  Ahm|,  using  the  above  expression  is  yet 
unaanageable  on  account  of  the  factorials#  but  now  a  stable 
recursive  formula  is  available: 


nmk 


[  (  2n-4k+l )  (  n-2k-m+l )  ( n-2k-m+2  )  (  n-2k+m+l ) (n-2k+m+2 ^  f 

2k(2n-2k+l)[2n-4k+51* 

±  (3.36) 

*nm  ,k-l 


°o 

1  <  k  <  s 


-n  <  m  <  n 


n  >  0 


with 


X 


nmo 


1 


for  all  n  ,  m 


(3.37) 


The  above  recursion  is  easily  obtained  from  (3.35)  by 
expanding  the  factorials. 

Poraulas  (3. 34),  (3.35),  and  (3.37)  show  that  the 
ellipsoidal  haraonic  coefficient  vnm  equals  the  spherical 
coefficient  of  the  saae  degree  and  order  plus  a  linear 
coabination  of  spherical  coefficients  of  lower  degree  (and 
saae  order)  .  The  ratio  B/E  for  the  earth  is  approximately 
<r.ssl2,  hence  froa  equation  (3.35),  the  coefficients  , 

It)  0  in  this  coabination  are  generally  auch  less  than  1  for 
low  values  of  n .  However,  increasing  values  of  n  coapensate 
the  rapid  decrease  of  1/o’etH  thus  also  slowing  the  rate  of 
decrease  of  Xww.t<  with  k ;  and  in  fact,  they  are  generally 
not  aonotonic,  since  tor  n>  4 <ra%2,  X*0|  >  1.  Therefore,  the 
larger  degree  ellipsoidal  harmonics  nay  have  considerably 
more,  or  less,  power  than  their  spherical  counterparts, 
depending  on  how  the  spherical  coefficients  conbine  to  fora 
the  ellipsoidal  coefficients.  Equation  (3.34)  also  shows 
that  a  finite  number  of  spherical  coefficients  generates  an 
infinity  of  ellipsoidal  coefficients;  thus,  if  the  function 
is  band  limited  in  spherical  coordinates,  its  ellipsoidal 
spectrua  is  infinite.  Bowever,  in  this  case  all  ellipsoidal 
harmonics  of  degree  higher  than  the  highest  spherical  degree 
are  linearly  dependent  on  the  lower  degree  harmonics. 

Because  the  equatorial  and  rotational  symmetries  are 
retained  when  transforaing  to  ellipsoidal  coordinates,  a 
spherical  series  of  even  and/or  odd  zonal  harmonics 
transforms  into  an  ellipsoidal  series  of  even  and/or  odd 
zonals,  respectively;  this  is  also  obvious  froa  equation 


(3.34)  . 

The  transf o nation  froa  ellipsoidal  to  spherical 
hareonic  coefficients  is  derived  similarly.  Substituting 
(3.26)  into  the  density  integral  (2.31); 


TCsrm(i>"  ^r-  t<2n+1> 

®  2n“2k(n-k) !  ,,,  .  .  N  «  ,*  u  ^ 

k!0<2n-2£+iynn  ///  n(u*6»x)  V2k,|m|  (1ir> 


•  V2k,m<5^>  * 


(3.38) 


which,  with  (3.29),  becoees 


(n+m) ! (n-m) !  (  E  .n+1  ?  (-l)k(n-k) ! (2n-4k+l) ! i 

2n+!  k£0  k!  (n-2k) !  (2n-2k+l) ! 


2n-4k+l _  1  - 

(n-2k-m) ! (n-2k+m) !  un-2k,m 


-  ?  L  G  s=  rp~lml  ] 

,1.  Lnmk  un-2k,m  »  s  1  2  J 


(3.39) 


where 


(-1)  (n-k)!(2n-ik+l)! 
k! (n-2k) ! (2n-2k+l)T 


C2n-4k+l)(n-m) ! (n+m) !  & 

f2n+l)  (n-2k-m) !  (n-2k+m)  r 


(3.40) 


or,  recursively, 


L  ■  1  ,  for  all  n,m 
□mo 
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3  (2n-2k+3) _ 

'nmk  2k(2n-4k+3)(2n-4k+5) 


(3.41) 


,  ( 2n-4k+l ) ( n-2k-m+l ) ( n-2k-m+2 ) (n-2k+m+l ) ( n-2k+m+2 

[ - 2n-4'k+5 - 


‘  Lnm,k-1 

1  <  k  <  s  ,  -n  <  m  <  n  ,  n  >  0 

Comparing  this  with  the  transformation  from  spherical  to 
ellipsoidal  coefficients,  the  general  coaaents  aade  for  the 
latter  obviously  apply  here  as  well.  It  is  equally  obvious 
that  each  transforaation  is  the  inverse  of  the  other  since 
the  set  of  harmonic  coefficients  is  unique  for  a  given 
function. 

Hith  the  renormalization  of  Q,,m  as  in  (3.32),  these 
functions  become  computational  tractable.  Proa  (Hobson, 
1965,  p.108)  and  with  the  usual  normalization  (1.14), 


(2n+l)(n-m) A 
Sn(n+m) !  J 


2nn!(n+m) 5 
(2n+lV! 


(y  2  -1)‘ 


im 


1_ 

y n+m+1 *  (3.42) 


F( 


n+m+2 


n+m+1 
~2 — 


n+- 


where  is  any  complex  number  with  \/*.\  >  1;  and  F  denotes 
the  hypergeoaetric  series.  Using  the  explicit  expansion  of 
P  (Abroaowitz  and  Stegun,  1970,  p.556)  ,  /a.  =io>,  and  equation 
(3.52),  we  find  for  the  function  Snm 


S  (o  p )  =  (1  +  -4 
nm  P  o 


im 

>  <• 


■)n+1  [1 


(m+n+l)(m+n+2) 
2  1 !  ( 2n+3 )  * 


1  (m+n+l)(m+n+2)(m+n+3)(ra+n+4)  l 

$  '  1  1  i|  “  •  •  ♦  J 

°  P  21-  2 !  < 2n+3)  ( 2n+5)  0  P 


He  note  that  as  E ■+  0  (ellipsoidal  coordinate  system 
degenerates  into  the  spherical  coordinate  system),  up-»rf, 
°V  oo  ,  P  -*  1,  o;/ov*a/u(,  -+  8/r^  and  hence 
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S„m  (crp)  -*  (B/r^)"4'  .  Also  -*  0  for  k  >  0,  so  that 

-•  vMW<  and  the  ellipsoidal  series  (3.33)  reverts  to  the 

spherical  series  (2.30)  ,  as  it  must  for  E=0 . 

For  later  use,  the  functions  are  similarly 
re normalized.  Let 


Snm<a)  “  (-i)”  /em(2n+i)(n-m)!(n+m)!  - -IJli p  (io)(3.«4) 

o?  in(2n+l)!  nm 


then  the  reciprocal  distance*  in  terms  of  Snm  and  is 


T  “  T  Jo  J.n  ^  S"I”>I<<’P)  5n|m|(0)  ^nm^p’^p^* 

•?__(«,  X) 
nm 


Bith  (Hobson*  1965,  p.95) *  we  find 


«  ,_■*  _  ,,  1  >>iin  ,  a  Nn  ^  (n-m)(n-ra-l)  1 

finm(a)  "  (1  +  ^}  il  +  —Sm=T) -  <3*46) 


which  is  a  finite  sum,  the  last  term  being  0  (<y”<"*m1)  . 
the  recursion  formula  for  the  Legendre  functions: 

P__(io)  ■  lo  o.  i  .  i  ( io)  -  B_  „  _  P_  0  _( ia) , 
nm  n-i,m  n-i,m  n-2,m  n-2,m  1 

n  >  2  ,  0  <  m  <  n  -  2 


Using 


5  i<i«)  “  io/2n+T  P  «.  n  idcr)  ,  n  2  1 

n,n-i  n-ifn-i 


(3.47) 


Pnn(io)  -  -IATWtF  Vl,n-l<10>  •  "  2  2 


where 
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a 


n-l,m 


./(2n-l)(2n+T7 

v(n-m)(n+m> 


0 


n-2,m 


_./(2n+l)(n+m-l)(n-m-f7 
V  (2n-3)  (!n+m)(n-ni) 


(3.48) 


we  obtain  a  stable  and  aacnrate  recursion  for  the  KHM : 
t?  /  — \_^s  /  —  \  i  ( n  +m-  l)(n-m-l)  1  s  /  _» 

Rnm  °  F;  Rn-l,m(a>  +  2^HT2n-T5  F\  Rn-2,m(a) 


n  >  2 


0  £  m  $  n  -  2 


R_  n  .(o)  =  Rn  -  n  .(a)  ,  n  * 
n,n— 1  Oo  n-l,n-l 


Rn  „<<» 

n  ,n 


S„  i  „  .(a)  ,  nil 
„2  n-l,n-l 


(3.49) 


3.2  Th£  Ellipsoidal  Series  o£  th$  gravity  Mi  Beiqht 

Anoaalies 

The  aoove  fora alas  establish  the  transforaation  froa 
spherical  to  ellipsoidal  haraonic  series  (and  vice  versa) 
for  an  arbitrary  (Newtonian)  potential  in  the  regions  where 
the  series  converge  tinifornly.  Utilizing  the  powerful 
theorea  that  haraonic  functions  and  Newtonian  potentials  are 
equivalent  (Kellogg,  1953,  p.218) ,  the  transforaation  of 
series  coefficients  of  the  gravity  anoaaly  and  height 
anoaaly  follow  ianediately  if  we  retain  their  definitions 
based  on  the  spherical  approziaation  (equations  (2.42)  and 
(2.57)).  Because 


r  Ag  *=  -r  -~-B  -  2T  » 
p  °Bp  p  3rp  ‘“p  * 


1 

n 


7Bt  p 


(3.50) 


are  both  haraonic  functions  in  the  exterior  space, 
application  of  equations  (3.33)  and  (3.34)  to  equations 
(2.50)  and  (2.58)  yields 


18<VVV  ■  ?rj0  J.n5»H<V  * 


nm  ?nm<VV  <3-5’> 
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with 


0 

Tnm  '  Jo  X“k 


.  s  -  [2* 


(3.52) 


and 


2  oo  n 


CP  nl  J-nS|1llll|(Cl|)  ?nm  ?ni"(ViP)  (3.53) 


where 


^nm  =  k|Q  Xnmk  zn-2k,ra  »  s  “  [2?3  <3-54> 


and  where  SHWI(  <rp  )  and  >ww,k  are  given  by  (3.43)  and  (3.35)  . 

The  coefficients  ,  or  I**.  ,  do  not  represent  the 
ellipsoidal  spectraa  (as  defined  by  (3.11))  of  the  gravity 
anoaaly  or  height  anoaaly  since  the  above  suas  are 
prenultiplied  by  r„  ,  which  is  a  function  of  u*,4P  and  hence 
not  constant  on  any  ellipsoid.  Finally  we  note  that  y.m  =0 
for  w=-1 ,0,1,  since  gln,s0;  also  if  g„  =0  and  z  ..  =0 ,  then  the 
corresponding  zero  degree  ellipsoidal  coefficients  vanish  as 
well. 


3.3  The  Derivation  of  the  Error  Series 

Using  the  fornnlation  of  the  disturbing  potential  as  an 
integral  of  a  generalized  density  layer,  the  downward 
continuation  error  of  the  ellipsoidal  harnonic  series  is 
similarly  derived  for  the  gravity  anoaaly.  Atteapting  a 
corresponding  developnent  for  the  potential  requires  the 
expansion  of  the  kernel  E,  above  and  below  the  bounding 
ellipsoid,  into  ellipsoidal  haraonics.  A  starting  point  for 
such  a  derivation  night  be  the  spherical  series  (2.115),  but 
this  aay  prove  to  be  a  fornidable  task.  It  is  not  pursued 
here  because  the  series  divergence  has  the  greater  relative 
effect  on  the  series  for  the  gravity  anoaaly. 
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Thus  consider  the  equation 


Hr 


//  ,X)  K  dm 

CD 


(3.55) 


where  d<o-sini  dS  dX  (ellipsoidal  coordinates)  and  R  (the 
notation  is  aodified  to  avoid  the  conflict  with  the  linear 
eccentricity)  is  the  saae  kernel  as  before  (equation 
(2.85)),  bat  in  ellipsoidal  coordinates;  that  is,  only  v 
has  changed  to  include  the  Jacobian  of  the  spherical  to 
ellipsoidal  coordinate  transformation.  Using  the  sane 
definition  of  the  gravity  anonaly  (equation  (2.42)),  we  get, 
as  before  (see  (2.91)), 


^Ep  ~  4ir 


=-  / /  v<6  iX  )  rgT^>dcD  (3.56) 

n  /a  *  n 


p  at 


By  the  addition  theorea  (1.9) 


-  i  4*  !  Sl,«p.>p)  13.57) 

p  p  m=-l  H  * 


Since  Tlw,  ( 8P,  Xp)/r£  satisfies  Laplace's  equation  and  is 
regular  at  infinity,  it  is  a  potential;  therefore,  the 
transf or nation  equations  (3.33) ,  (3 .34)  fron  spherical  to 
ellipsoidal  harmonic  series  apply: 


Ylm(0P 


eo  n 

Jo  JLJ-mV  n‘”‘  (3-58) 


wh  ere 


^nt 


s 

k*0  *ntk  ^n-2k,t 


(3.59) 


and  the  qwt  are  the  corresponding  spherical  haraonic 
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coefficients : 


qnt  = 


0 

R2 


if  n  *  1 

if  n  =  1  and  t 


-1,  0,  1 


(3.60) 


Inserting  (3.60)  and  (3.59)  into  (3.58),  we  find 


■  X  Sn|»|(«>  Snm  fnm<VV  (3-61> 

n  odd 

Finally,  by  substituting  (3.61)  into  (3.57)  and  considering 
the  transf or nation  between  inner  harnonic  functions, 
equation  (3.26) ,  we  obtain  the  foraula 


n 


rs  cost(i  _  1 

“ -  ™nL0  mi_nu““  "M-IV”  °n I m 1 1 1 » | m | 


rJ -  Tm  I  I  Urn*>  lmf(6»x)* 

P 


(3.62) 


•  Y  <  6  ,  X) 
nm  p  p 


where  tne  explicit  expression  for  uHm  is  not  difficult  to 
find,  but  of  no  consequence  in  the  present  derivation;  it  is 
noted,  however,  that  u,m*1,  for  all  n.  Combining  (3.62) 
with  the  ellipsoidal  series  expansion  for  the  reciprocal 
distance,  equation  (3.46),  and  inserting  this  into  (3.56) 
results  in 


it!p  *  IWT  &  Sn|m|(V 

*  P  n®0  m*-n  1  1  *  w  _  (3. 


63) 


•  Y 


nm'  p'  p 


where  the  point  P  is  located  above  the  bounding  ellipsoid 
(u,>  b)  and 
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“am  "  T[S^I  Sn|m|<0)  ?mn(S  -X >  *  Tunm  Bl,l»l<9>'  <3-6,> 

(Mote  that  =0.)  It  the  earth*s  surface  is  an  ellipsoid 

(with  seaiainor  axis  b ,  e.g .  b=6 356755 .288  a  for  the  aean 
earth  ellipsoid)  on  which  the  gravity  anoaaly  is  a  known 
function , 


Agp  =  rp  J[0  mJ_n  §n|m|  (T5  ^nm  ?nm(6p*Xp)  <3*65> 
n#l 

then  in  this  case  the  integral  equation  (3.63)  can  be 
solved,  since  a  coaparison  with  (3.65)  shows  that 

*nm  *  WH  V(5-X>  «,X)du 

03 

=  Wfj  v(5’X)  <25+1  ^nlml  055  5nm(5  ’X  }  S1 J  m[  (°  > 

^lm(6,X)lda>  (3.66) 

e  E&l  *n\m(3'>  §n|m|(5)  Dnm "  Tunm  Rl,|m|(5)  Dlm 


where  O’  =b/E  and 


v(B,S,X)  =  |  J  S  |  |  (5)  Dnm  ?nm(s.x)  (3.67) 

n*»u  m=-n  ' 


Since  u,m=1  and  y,*,=0,  D,m  cannot  be  deterained  froa 
equation  (3.66);  we  aay  assuae  D0o  -Dtfr  =0 .  Then 


2n+l 


nu> 


Rn|m|  <iT>  sn|m|<a)  Y”" 


(3.68) 


For  points  below  the  bounding  ellipsoid,  consider  the 
regions 


Ill 

“i  “  {(5, A)  /  up  >  b}  ,  o)2  =  { (6  ,  X  >  /  Up  <  b}  (3.69) 

If  (  ,  the  reciprocal  distance  is  expanded  as 


n 


l“  =  T  Jo  §nM«”fWV  ?n„,<5-X>' 


ynm(6n,A^) 

run  p  p 


(3.70) 


Substituting  this  and  the  series  for  r,  cos y /rl  (3.62)  into 
(3.56)  yields  p 


+ 


R 

4nr 


00 


-  I 

p  n=0 


1 


m=-n 


//  v(6,A) 

0)2 


l2n+I  Sn|m|(op>  §n[m|(o>‘ 


•Yn.n(6«X> 


(3.71) 


-  -Y-  Sn|m|<V  Hl,|m|<0>  »»)*•  *„m<VV 

Mote  that  ( r  =u/E  is  a  function  of  (  6  ,  A  )  ,  since  u  is  the 
coordinate  of  a  point  on  the  earth's  surface. 

The  downward  continuation  error  of  the  ellipsoidal 
harnonic  series  for  &  g  is  then  the  difference  between 
(3.63),  truncated  at  n=n,  and  the  true  series  (3.71).  After 
several  simplifications,  in  which  (3.64)  is  duly  considered, 
this  error  can  be  expressed  as 


n  n 

e(Ag  )  »  l  l  d _ ( o  )  Y__(S  , A  )  + 

p  n=0  m--n  nm  p  nra  P  P 


+  l  l  d*  (a)  ?  (6  .A  ) 

n-5+1  m— n  nm  p  nm  p  p 


(3.72) 


* 


where 
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■WV  “  }J  V<S,X)  tBB|B|<op)  E„|m|<°>  (3.73) 


-  Sn!m(VSnM(<,)1?nn,<^>d“ 


"nn/V  -  <WV  '  ^  Sn| ml  <°p>  \m 


(3.7*) 


being  the  ellipsoidal  coefficients  of  the  outer  series 
of  r,Ag,.  In  view  of  (3.43)  and  (3.46)  the  kernel  of 
(3.73)  can  be  expanded  as  a  series  in 


i  «  1  -  2a 

a 


(3.75) 


By  using  the  series  expressions  for  Sn>M(<r)  and  B^tr), 
naaely  equations  (3.43)  and  (3.46)  ,  and  collecting  like 
powers  of  w  and  arp**,  it  is  found  (through  a  lengthy  and 
tedious  derivation)  that 

Snn1<o  )  R  <o)-S  <o)  R__(o_)  - 

nm  p  run  nm  nm  p 

=  (2n*J)^lr5  (  1+  Tii-i-  till-;  +  0(i  ))  + 

°P  °P  dP  °P 


where 


+  W2(  T2lJL  +T22-L  +  0(i  ))  + 


+  W*Cr$«+Tjl-~  +  T22-V  +  0  ( — t  )  )  + 
°P  °P  °P 

+  +  0(~.  ))  + 

°P  °P 

♦  0(w* ) ] 


Tn  -  -1  ,  1x2  -  1  ,  Tji  ■  1  ,  T22  -  "2 


(3-76) 


.  n(n+l)  _  __  mV(2n+l)(n+l) 

T»#  ■  T*#  ■  — -  ,  t|j" - — - 

t 1 2  ■  2  +  i(na+n+m*)  ,  t*2  ■  1-iu 


(3.77) 
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As  with  the  total  downward  continuation  error  in  spherical 
haraonic  series,  the  linear,  as  well  as  qaadratic,  tens 
teras  in  w  son  to  zero  for  points  on  or  above  the  earth's 
surface.  The  expression  (3.76)  provides  two-  to  three-digit 
accuracy  for  n  =  300,  *=0,  150,300  and  o>  =  12. 180886 , 
v  =12.193239  which  represent  the  worst  situations;  that  is, 
when  the  point  of  computation,  P,  and  the  point  of 
integration  are  farthest  apart  in  tens  of  the  coordinate  a 
(then  the  ratio  «v/cr  is  least  and  5  is  largest)  . 


3.4  The  Bespits  of  the  Haaerical  Analysis 

Since  the  infinite  ellipsoidal  haraonic  series  of  the 
gravity  (or  potential)  converges  with  certainty  only  outside 
the  bounding  ellipsoid,  the  evaluation  of  the  truncated 
series  at  the  earth's  surface  is  associated  with  a  downward 
continuation  error.  However,  we  should  expect  the  error  to 
be  saaller  than  for  the  spherical  series  since  the 
penetration  into  the  region  of  (probable)  divergence  is 
generally  not  as  deep,  especially  in  the  polar  areas. 
Although  this  is  alaost  obvious,  it  can  be  verified  by 
exaaining  the  divergence  of  the  zonal  series  corresponding 
to  the  siaple  density  distributions  of  section  2.1.  For 
exaaple,  Pigure  19  shows  the  differences  between  the  partial 
suas  of  the  ellipsoidal  series  and  the  true  value  of  the 
gravity  anoaaly  evaluated  at  the  point  (r* =6357200  a, 
®,=7£S).  The  eccentricity  of  the  coordinate  systea  was 
taken  as  2=450000  a  so  that  the  evaluation  point  lies  below 
the  bounding  ellipsoid:  -5300  a  for  the  equatorial  disk  and 
-3500  a  for  the  serrated  ellipsoid.  The  ellipsoidal 
haraonic  coefficients  were  deterained  by  applying  the 
transforaatian  (3.52).  Clearly,  in  contrast  to  Pig.  5,  the 
effect  of  divergence  is  sore  subdued,  becoaing  noticeable 
only  when  n  >  1800. 

The  above  expectations  are  not  realized  when  coaparing 
the  truncated  ellipsoidal  series  of  the  gravity  field 
generated  by  the  density  layer  of  section  2.3  against  the 
corresponding  truncated  inner  spherical  series.  The 
differences  between  the  partial  suas  for  n  =300  can  be  orders 
of  aagnitude  larger  than  the  values  listed  in  Tables  3  and 
4.  The  only  adaissible  conclusion,  that  this  is  a 
coaparison  of  incoapatible  spectra  and  is  therefore 
aeaningless,  reemphasizes  the  inherent  danger  in  the 
coaparison  of  partial  suas  of  different  series  representing 
the  saae  function. 
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Figure  1*»:  Partial  suas  of  ellipsoidal  haraonic  series  of 
a  g  tor  eguatorial  disk  and  serrated  ellipsoid 
density  distriJbntions  (Figures  i  and  3)  ainus 
corresponding  true  values  evaluated  at 
rp  =6357200  m,  0f*7!5. 
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For  the  numerical  study  of  the  downward  continuation 
error,  we  can  specify  the  density  nodel  by  (3.68)  in 
conjunction  with  the  coefficients  (2.126).  The  quantities 
were  set  equal  to  g„„ ,  recognizing  that  a 
different  gravity  nodel  ensues  and  that  the  correlation  to 
the_ topography  is  no  longer  given  by  (2.131).  This  choice 
of  is  equivalent  to  the  assuaption  that  the  function 
iP h  g*  (approx iaa tel y  on  an  ellipsoid)  was  analyzed  in  the 
ellipsoidal  coordinate  systen  (  <»>  )  to  yield  Hg«...  (The 
degree  variances  of  D**,  (3.68)  and  of  (2n  +  1)gnw,,  equation 
(2.94),  differed  only  in  the  third  digit.)  Horeover,  to 
avoid  unnecessary  complications,  the  sane  surface  aodel  was 
adopted,  but  the  grid  coordinates  were  identified  as 
ellipsoidal  coordinates:  S{=  in6  ,  where  aS=  0!6.  Hence 

usl  *  Vrs!i  -  Si,f6Ui  (3.78) 


The  seaiainur  axis  of  the  bounding  ellipsoid  (with 
E~521854 .4492  a,  see  (3.9))  was  found  to  be  uw  =6 363096. 071 
a.  The  surface  deviates  froa  this  ellipsoid  by  an  BUS  value 
of  6114  a ;  the  deviation  is  around  6000  a  whether  at  the 
pole  or  equator.  Other  choices  of  the  linear  eccentricity  E 
aay  produce  a  closer  overall  fit  to  the  surface;  for 
exaaple,  with  £=523836.8873  a  (a=6378140  a,  f*1/296),  the 
bounding  ellipsoid  is  closer  to  the  surface  at  the  poles, 
bat  aore  distant  at  the  equator,  with  a  total  BBS  deviation 
of  6085  a. 

The  evaluation  of  the  coefficients  d„m  of  the  downward 
continuation  error  of  ellipsoidal  series  is  alaost 
identically  perforaed  as  for  the  spherical  series.  The 
expansion  (3.76)  ainus  the  linear  and  quadratic  terms  in  w 
is  substituted  into  (3.73),  which,  in  turn,  is  discretized 
according  to  the  assuaption  that  the  surface  aodel  consists 
of  ellipsoidal  ooapartaents  delineated  by  the  coordinate 
lines  &  ^constant,  X=constant,  and  that  v  is  a  step 
function  constant  within  each  coapartaent. 

in  contrast  to  the  determination  of  the  error  in  the 
spherical  haraonic  series,  one  value  of  up  (u„=6356800  a) 
sufficed  to  yield  an  ellipsoid  with  an  adequate  supply  of 
points  within  100  a  above  the  surface  aodel,  in  all  regions. 
The  BHS  values  of  the  downward  continuation  error  series, 
truncated  at  5=300,  at  these  points  are  shown  in  Table  8  for 
the  regions  of  Table  1  (where  the  latitudinal  ranges  now 
refer  to  the  reduced  latitude) .  Also  shown  are  the  numbers 
of  total  points  on  which  the  HNS  values  are  based,  as  well 
as  the  aaxiaua  absolute  value  in  each  group. 
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Table  6:  Hhs  downward  continuation  error  of  ellipsoidal 

haraonic  series  of  gravity  anoaaly  in  regions  of 
Table  1  (latitude  ranges  refer  to  reduced  latitude) 
at  points  on  the  ellipsoid  Up=6356800  a 
(£=521854.4492  a)  and  above  the  surface  aodel  no 
aore  than  100  a.  flaxiaua  absolute  values  for  each 
region  are  given  parenthetically. 


Beg ion 

has  £(ag)  (aax  £(ng))  ,  ,*gal, 

A:  o!3  -  359#.7 

0  i  n  i  300 

No.  of 
points 

M 

No.  of 
points 

S 

I 

74  10 

1.9x10"'’  (3.2x10  *  ) 

83  55 

2.7x10'"  (4 .  lx  10' J  ) 

11 

500  3 

1.7x10'''’  (1.  1x10‘*  ) 

9328 

.94x10*"  (.45x10  *  ) 

111 

7282 

1.1x10  "  (.96x10**  ) 

20  71 

1.3x10  "  (.63x10**  ) 
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The  nuaerical  results  are  Halted  to  the  demonstration 
that  the  effect  of  divergence  (h=300)  of  ellipsoidal 
haraonic  series  near  the  earth's  surface  is  considerably 
less  pronounced  than  in  the  case  of  series  of  spherical 
haraonics.  Mote  also  the  essential  uniforaity  of  the  errors 
over  all  the  latitudinal  ranges. 


4.  Corrections  to  £he  Spherical  Approxiaations 


The  ellipsoidal  correction  expounded  by  Lelgeaann  (1970) 
and  Horitz  (1980)  aaend  the  spherical  approximation  by 
accounting  for  the  general  ellipticity  of  the  earth's  shape. 
The  spherical  approximation  of  the  relevant  geodetic 
quantities,  ng,  if  ,  defined  by  (2.42)  and  (2.57),  however, 
do  not  conform  precisely  to  Moritz's  definitions.  To 
achieve  the  spatial  spherical  harmonic  expansion  of  the 
height  anoaaly  (equation  (2.56)),  as  well  as  the 
corresponding  ellipsoidal  expansion,  the  noraal  gravity  was 
equated  with  the  gravity  produced  by  a  homogeneous  ball  of 
aass  H,  instead  of  the  conventional  average  value  of  y  over 
the  spheropotential  surface  (the  latter  being  assuaed  by 
Horitz) .  Furthermore,  the  spherical  approximation  of  the 
gravity  anoaaly,  according  to  Horitz  (ibid.,  p.425),  is 
defined  by 


Tp  -  J!  Ag°  S(*)  do  (4.D 

*  rt 


which  represents  the  solution  to  the  (third)  boundary  value 
problem  if  the  bounding  surface  is  a  sphere;  S(y-)  is 
Stokes*  function.  In  equation  (4.1)  T  is  the  actual 
disturbing  potential  with  no  spherical  approxiaation,  and 
both  T  and  ag°  are  functions  on  the  ellipsoid  that 
approximates  the  earth's  surface,  with  seaiainor  axis,  say 
b.  Using  coordinates  for  which  the  ellipsoid  is  a 
coordinate  surface,  for  exaaple,  ellipsoidal  coordinates 
(u,  i  ,  X  )  (Horitz  uses  geodetic  coordinates),  the  angle  -y 
loses  its  usual  geoaetric  meaning  since  it  is  defined  by 


cosp  *  cos  6  cos  6p  +  sind  sindp  cos(X-  Xp)  (4.2) 

where  (  £,  X  )  and  ({,  ,X,)  are  points  on  the  ellipsoid .  y 
aay  be  interpreted  as  the  central  angle  between  the 
projections  of  ellipsoidal  points  onto  the  unit  sphere 
according  to  the  correspondence  (3.12) .  Expressing  T  and  S 
in  teras  of  harmonic  functions. 
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T(b*VV  *  J  l  V  ?nn,(VV 

*  n=»0  m=-n  v  v 


(4.3) 


00  Q 


SW,X,«p,»p.)  -  nI2  J_n  5^  yi,»  *nm(VV 


directly  yields 


oo  n 


n-1 


Ag"(b,{,X)  -II  ^  An„?nm<6,X) 
ns2  m«-n 


(4.5) 


Equation  (4.5)  obviously  differs  fro*  (3.51)  because  of  the 
latter's  dependence  on  r .  Bence  the  ellipsoidal  corrections 
of  Boritz  (1980 ,  pp. 318, 327)  to  ?  and  ag  are  not 
applicable  here. 

When  abandoning  the  spherical  approximations  (2.42)  and 
(2.57),  the  siaple  spectral  relationship  between  the 
potential  and  gravity,  or  height,  anoaaly  is  lost, 
nevertheless,  to  the  approxiaation  developed  below  and 
knowing  their  relationship  in  the  space  doaain,  the  latter 
are  still  representable  as  series  involving  the  spherical 
harmonic  functions. 

In  the  following,  we  will  build  on  the  premise  that  the 
disturbing  potential  is  known  to  any  desired  accuracy,  for 
example,  as  a  series  of  spherical  or  ellipsoidal  haraonics. 
The  corrections  to  the  spherical  approxiaation  of  the  height 
anoaaly  and  gravity  anoaaly  will  be  derived  to  an  accuracy 
determined  by  the  neglect  of  teras  involving  the  fourth 
power  of  the  first  eccentricity.  Furtheraore,  derivatives 
along  spheropotential  surface  normals  and  along  ellipsoidal 
noraals  are  not  distinguished  (they  are  identical  on  the 
reference  ellipsoid,  if  it  is  an  eguipotential  surface  in 
the  normal  gravity  field).  Thus  the  height  anoaaly  is 


C 


P 


~[T  + 

P 


2  32U 
CP  3h2  Iq 


+ 


•  3 


(4.6) 


where  ‘V^h  is  the  directional  derivative  along  the 
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ellipsoid  noraal.  The  vertical  gradient  of  the  noraal 
gravity,  ^  y  /  ^  h*  "Su/ *2  h1- ,  is  on  the  order  of 
2  y /r » 3x10*‘ s*1  ,  so  that,  with  <  100  a,  its  oaission 

causes  at  aost  an  error  of  3  aa.  In  order  to  achieve  0.5  ca 
accuracy  in  the  height  anoaaly,  the  noraal  gravity  aust  be 
accurate  to  about  50  agal,  so  that  in  view  of  the  above 
gradient  of  about  0.3  agal/h,  Vv  aay  be  substituted  for 
in  (4.6)  .  Similarly,  the  gravity  anoaaly  in  its  aost 
rigorous  fora  is 


C2  + 


(4.7) 


where  W  is  the  earth's  gravity  potential,  0  is  the  noraal 
potential,  and  Jd/3u  is  the  derivative  along  the  pluab  line, 
i.e.  the  gradient.  Terns  of  second  and  higher  order  can  be 
neglected,  causing  an  error  of  at  aost  1.5x10~3  agal,  since 
'3vy/‘3h'ss  6  k  h/b1*  ~  1.5x10'"’  a*'  s‘l  .  Purtheraore,  as  this 
also  gives  the  change  in  the  vertical  gradient  of  y  with 
height,  the  valae  of  Ag  changes  by  no  nore  than  1.5x10_1 
agal  if  ly/3h  is  evaluated  at  P  instead  of  Q. 

Pig.  20  shows  the  direction  of  the  plnab  line  with 
respect  to  the  orthogonal  directions  of  dh,  Zd<p  ,  Ncos?  dl 
at  a  point  P;  Z  is  the  aeridional  radius  of  curvature  of  the 
ellipsoid,  N  is  the  radius  of  curvature  in  the  priae 
vertical,  and  (p  is  the  geodetic  latitude.  If  0  is  the 
total  deflection  of  the  vertical  (angle  between  dh  and  dH) 
with  coaponents  J  and  »j  ,  then,  since  these  are  snail 
angles,  the  direction  cosines  of  dB  with  respect  to  the 
noraal  directions  are  cos  9 ,  S  #  and  >|  .  Bence 


TR  m  cosQ  +  *  zV<h  +  YcoSlax 


(4.8) 


Considering  Pig.  21,  the  directional  derivative  along 
the  ellipsoid  noraal  is 


■  oos*  -£r  '  si“'"  T5e 


(4.9) 
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where  y  is  the  angle  between  dh  and  dr.  The  ellipsoid 
nornal  intersects  the  egoator  at  a  distance  s  from  the 
center;  we  have 


aeVl-e2  sln0 
/I  -  e2  sin20 


a  /l  -  ez 
r  =  - - - 

/l  -  e2  sin2 9 


(4.10) 


Using  the  law  ot  sines  on  the  triangle  OBQ,  we  obtain 


£. 

s 


=  cottp  cos0  +  sin0 


(“.11) 


Hence  with  (4.10) 


tan  ip  - 


e  2  sinOcosO 
1  -  e2sin20 


(4.12) 


Pythagoras*  theorea  then  easily  furnishes  the  expressions 


sin^  =  e2sin0  cosfl-  ie4  (e2- 2)  sin30  cos6  +  0(e6  )  (4.13) 


cosip  =  1  -  ie4sin20  +  2*.  e4  sin40  +  0(ee>  (4.14) 


Hence,  neglecting  terns  of  0(e**)f 


TB’-y  -e2sin9  cosQ  FT0  (4. 15) 


Note  that  e4  depends  on  the  coordinate  surface  u=b  under 
consideration,  but  to  the  accuracies  involved  here,  it  can 
be  treated  as  a  constant  for  points  near  the  earth's 
surface. 


The  nornal  potential,  a,  is  given  as  a  series  of 
spherical  haraonics  by  Heiskanen  and  Boritz  (1967, p. 73): 
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U  -  ~[1-  l  J2n(T)2n  p2n(cos0)1  +  ^r2sirfe  (4-16) 
n=l 


where  Jt..*=o  (e2” )  and  oj  is  the  rotational  speed  of  the 
earth:  to1  »  5x10"1s*t  .  The  nornal  gravity,  being  the 
gradient  of  0,  is  therefore 


Y  =  _  JE  =  _2E-  +  e2sine  cos  0~ 
y - 9h  §r  e  sme  cos  er30 


3-^J2(-S-)2  P2(cos9)  -  u»2r  sin*0(l  -  e*cos?0) 


+  0(e%  ) 


(4.17) 


He  note  that  co*r »  0 .03m/s"ts  0.003 y  »  0(ev),  so  that  the 
terms  with  cox  rev  can  also  be  neglected.  From  (4.17),  the 
normal  derivative  of  y  is  found  to  be 


-  e2sin0  cos0- 


=  -  ~tY  -  3-p- J2(^-)2P2(cos0)  + -|-u)2r  sin20]  +0(e' 


(4.18) 


Hence 


1  3y 
y  Its 


-  ~  +  (^-  J2(-fr)2P2  -au^irfP  +  0(ev)) 


(1  -  3J2(~)2P2  -w2r  sin^*  OCe*))' 


(4.19) 


.  _  ^-[1- 3J2(-p)2P2(cos0)  + -|-^,sin20] +0(e") 


Substituting  (4.17)  into  (4.6),  the  height  anomaly  becomes 
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Cp  -  ^^[l  +  3J2(^)2P2(cos9p>  +  ^sin%p]  +  0(e“  )  (4.20) 


The  correction  above  teaches  a  aaxiana  value  (0.33%)  at  the 
poles,  8,=0®  ,  180®  . 

Taking  note  of  (4.15)  and  (4.8),  the  gravity  anoaaly 
(4.7),  upon  substituting  (4.19),  becoaes 


4gP=-l^-e,sl%cos^^f5;+i  &ip  Vep+0(e>) 


-  -  H5  -  J-  T  -  e2sin0  cose  + 

3rn  rD  P  P  P  rn“°D 


P°°P 

+  (6J2^-P2(cose  )--|fPsin20  )  T  +e  +0(e*) 

p  v  ' 


(4.21) 


where  Tp=«p-U^  is  the  disturbing  potential  ana 

£P  (1  “  cosQp)  3H^  ~  ^p  Zpd$p  "  np  Np  cos^Xp 


=  (1  -  cosO_) 


P 
3  Wr 


,  i*B  -  E  3T°  -  „  \  TP  ,r  i02 

P  9hp  P  zp3*p  P  Np  c°P^9Jp  Vzp9* 


(4.22) 
3Up 


in  which  we  used  the  fact  that  the  noraal  potential  does  not 
vary  in  longitude.  The  derivatives  of  T  in  (4.22)  are  the 
components  of  the  deflection  of  the  vertical  aultiplied  by 
the  noraal  gravity,  and  the  derivative  of  U  is 


3  0  _ 


-|H—  -  j2  +  w2r)  sin9  cos  6 

R36  R 


Finally, 


3W  _ 

TFT  " 


g 


(4.23) 


(4.24) 
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so  that  an  approximate  upper  bound 
is 

Up!  *  le2pgpl  +  |cp  V 

+  |  £p( J 2  +  u2  R)  sin9 

Using  gays  9.8  m/s1-  and  5  *  as  10",  ve  find  U,l  <  0.34 
■gal;  the  dominant  ter*  is  the  last  one,  hence  this  bound 
varies  linearly  with  the  magnitude  of  In  view  of  its 

relation  to  the  deflection  of  the  vertical,  the  horizontal 
derivative  of  T  has  typical  values  of  *  £  *50  mgal,  with  an 
upper  bound  (£<!•)  of  300  mgal.  Therefore  the  third  term 
in  (4.21)  has  a  value  of  about  0.2  mgal  (at  most  1  mgal), 
while  the  last  two  terms  in  (4.21)  rarely  amount  to  more 
than  0.1  mgal.  The  formulas  (4.20)  and  (4.21)  are  valid 
anywhere  on  the  earth's  surface  or  above  with  an  accuracy  of 
about  e^p  <  5  ■■  and  e^og  <  1x10*1  mgal,  respectively. 
Equation  (4.21)  generalizes  the  correction  derived  by 
flolodenskii  et  al.  (1962,  p.212)  for  anomalies  on  the 
reference  ellipsoid  and  is  quite  unlike  the  correction  of 
fioritz  because  of  the  different  definition  of  spherical 
approximation. 

The  equations  (4.20)  and  (4.21)  were  developed  in  order 
to  provide  the  means  for  the  precise  evaluation  of  the 
gravity  and  height  anomalies  using  the  spherical  (or 
ellipsoidal)  harmonic  series,  whose  formulation  depends  on 
the  adopted  spherical  approximation,  namely  (2.42)  and 
(2.57)  . 


for  the  magnitude  of  z? 


+  n2  Y 
1  P  P1 


COS0  I 


(4.25) 


5.  Summary.  Conclusion.  Becoaaendation 


The  expansion  of  the  earth's  gravitational  potential 
into  a  series  of  spherical  haraonic  functions  has  long  been 
used  to  describe  it  on  a  global  basis.  The  question  of  the 
validity  of  such  an  expansion  at  the  earth's  surface,  though 
propounded  t rom  the  outset,  has  been  addressed  firaly  only 
recently  and  then  primarily  froa  a  purely  theoretical 
standpoint  (e.g.  the  Hunge-Krarup  theorea) .  while  these 
represent  iaportant  advances,  a  definitive  answer  has  yet  to 
be,  or  aay  never  be,  found.  In  the  practical  situation,  the 
infinite  series  is  necessarily  truncated  at  soae  degree  n. 
What  effect  the  (possible,  or  even  probable)  divergence  of 
the  infinite  series  at  the  earth's  surface  has  on  its 
partial  suas  has  received  only  a  "first-generation1* 
analysis,  kith  the  present  study,  we  have  taken  a  second 
look  at  this  effect,  but  the  subsequent  conclusions  aust  be 
carefully  phrased  and  are  necessarily  lacking  in  nuaerical 
specificity. 

The  downward  continuation  error,  being  defined  here  in 
connection  with  the  effect  of  series  divergence,  is  a 
deterainistic,  or  systeaatic,  error;  it  has  no  stochastic 
properties  (one  could  argue  this  point  if  the  convergence 
surface  is  itself  a  stochastic  process)  .  Therefore,  its 
assessaent  is  forthcoaing  if  the  true  value  of  the  field 
which  the  series  represents  is  known,  such  as  in  the  case  of 
the  siaple  density  uistr ibutions  of  section  2.1.  These 
distributions  were  not  designed  to  simulate  the  earth's 
density,  but  their  diaensions  and  the  nuaerical 
investigations  were  selected  with  the  terrestrial  situation 
in  aind.  It  was  found  that  the  partial  suas  of  the 
spherical  haraonic  series  for  the  potential  (and  gravity 
anoaaly) ,  evaluated  below  the  surface  of  convergence,  do  not 
show  signs  of  divergence  until  the  truncation  degree  is 
relatively  large.  Also,  Figure  8  suggests  that  the 
divergence  problea  aay  affect  the  geopotential  series  only 
if  the  truncation  degree  is  300  or  greater.  For  the  earth, 
the  requireaent  of  knowing  the  true  value  of  the  potential 
field  at  the  surface  is  difficult  to  aeet,  even  if,  for  the 
purpose  of  an  ad  hoc  analysis,  aodels  are  introduced  to 
represent  the  earth's  surface  and  gravity  field. 
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harly  works  were  based  on  modeling  the  earth's 
lithospnere  by  a  voluaetric  density  distribation  of 
essentially  constant  value.  In  this  and  the  present  use  of 
a  surface  layer  density  distribution,  the  true  values  of  the 
modeled  field  can  only  be  estiaated,  again,  by  a  truncated 
(inner)  series  expansion.  It  was  found  that  low  degree 
comparisons  of  inner  and  outer  series  of  the  voluaetric 
density  aodel  give  unrealistically  large  values  of  the 
error.  This  result  is  attributable  principally  to  the 
inadequacy  of  the  choice  of  the  aodel.  Furthermore,  since 
this  aodel  itself  depends  only  aarginally  (through  a  crustal 
density  of  2.67  g/ca1,  etc.)  on  the  actual  characteristics 
of  the  earth's  potential  field,  it  was  abandoned  for  a 
surface  density  layer  aodel  which  could  be  defined  (although 
not  optimally)  so  as  to  yield  a  reasonable  representation  of 
the  true  potential  field.  This  aodel  does  not  represent 
exactly  the  earth's  field  at  resolutions  greater  than  200 
km ,  since  randoa  harmonic  coefficients  were  generated  to 
fill  in  the  detail  to  a  resolution  of  67  ka.  However,  the 
coefficients  were  forced  to  decay,  with  degree  n,  according 
to  a  degree  variance  aodel  characteristic  of  the  true 
gravity  field.  The  aain  disadvantage  of  this  density  aodel 
is  the  inability  to  compute  the  true  values  of  the  potential 
on  the  surface  to  any  desired  accuracy.  A  aodel  for  which 
this  is  no  problea  consists  of  a  sufficiently  large  nuaber 
of  point  masses  distributed  globally  just  below  the  earth's 
surface  (to  ensure  the  divergence  of  the  series)  ,  as  well  as 
deeper  within  (to  generate  long  wavelength  power).  The 
difficulty  with  this  aodel  would  be  the  numerical 
determination  of  the  aasses  for  a  representative  potential 
field. 

The  downward  continuation  errors  depicted  in  Tables  3 
through  7  are  completely  insignificant  with  respect  to 
anticipated  measurement  accuracies  of  1  agal  and  10  ca  in 
the  gravity  anoaaly  and  geoid  undulation,  respectively.  For 
the  anoaaly  the  sua  of  the  haraonics  of  the  error  up  to 
degree  300  was  found  to  be  0.3  -  0.5  yugal  (1  /<gal  = 

10°  mgal)  near  the  equator,  about  20  >*-gal  in  the 
aidlatitudes,  and  80.  -  90.  >4 gal  in  the  polar  regions. 
Similar  minute  values  were  obtained  for  the  first  300 
degrees  of  the  error  in  the  height  anonalies:  0.3  -  0.5  y  a 
(1>ab  =  10*‘  m)  in  the  low  latitudes,  about  70  /a.  a  in  the 
aidlatitudes,  and  approximately  h00  /*.  a  near  the  poles.  Of 
course,  these  numbers  do  not  give  the  entire  error  since 
they  exclude  the  contribution  fora  degrees  301  to  o o, 
representing  a  truncation  effect,  i.e.,  the  neglect  of  terns 
of  degree  greater  than  S  of  the  inner  series.  We  can  expect 
the  usual  degree  variance  models  to  provide  a  fair  estimate 
of  this  effect.  This  expectation  is  rather  intuitive,  based 
on  the  near  sphericity  of  the  earth's  surface,  the  results 
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of  section  2.1,  and  cm  the  smallness  of  the  values  of  Tables 
3  and  4,  but  not  founded  on  deductive  reasoning;  the  same 
expectation  say  prove  to  be  erroneous  for  the  voluaetric 
density  aodel.  Using  the  aodel  (2.128),  the  BUS  truncation 
error  (3-300)  has  values  (for  point  estimates)  of  about  30 
agal  and  36  cm,  respectively.  Por  the  Gaussian  smoothed 
fields  (954  of  the  300-th  degree  harmonic  is  filtered) ,  the 
first  300  degrees  of  the  gravity  anomaly  error  approximately 
sum  to  0.02/*gal  (equatorial  region),  1.0  >cgal 
(■id latitudes)  ,  and  4.0  /<gal  (near  the  poles)  ,  with 
respective  values  of  0.02>cm,  3.0  /cm,  and  20/cm  for  the 
height  anomaly.  The  BUS  truncation  effect  is  approximately 
0.3  mgal  and  0.7  cm,  respectively.  Therefore,  the 
estimation  of  point  or  mean  gravity  anomalies  and  geoid 
undulations  (height  anomalies)  using  the  outer  series 
expansion  to  degree  300  anywhere  on  the  earth's  surface  is 
practically  unaffected  by  the  divergence  of  the  total 
series. 

Throughout  this  exposition  emphasis  has  been  on  the 
dissimilarity  of  partial  sums  of  inner  and  outer  series. 

Por  a  constant  radius  r,  =8,  greater  than  the  bounding  sphere 
radius  H,  the  coefficients  of  the  partial  sum  of  the  outer 
spherical  harmonic  series  for  the  potential  represent  a 
portion  of  its  spectrum  on  the  sphere  of  radius  Bp.  Since 
the  spheres  of  radius  rf <  B  pass  through  the  earth's 
interior,  the  coefficients  of  the  inner  series  (constant  rp) 
cannot  represent  the  spectrum  of  the  exterior  potential. 
Indeed,  for  the  density  layer  model  it  is  difficult  to  give 
an  interpretation  to  these  coefficients,  which  in  any  case 
vary  as  the  point  P  moves  on  the  earth's  surface.  By 
accepting  the  conclusion  that  the  truncated  outer  sseries 
(n  £  300)  can  be  used  without  concern  for  divergence  anywhere 
on  the  earth's  surface,  we  also  cannot  identify  the  outer 
harmonic  components  as  spectral  constituents  of  the 
potential  (or  gravity  anomaly)  on  the  surface.  (Mote  that 
the  harmonics  should  be  evaluated  on  the  actual  surface  of 
the  earth,  and  not  on  some  mean  earth  sphere.)  Therefore, 
any  evaluation  of  the  outer  series  must  be  accompanied  by  an 
unambiguous  statement  regarding  the  guantity  being 
estimated.  Por  example,  by  introducing  the  Gaussian 
average,  we  attempt  to  eliminate,  or  filter,  the  high-degree 
information,  so  that  the  inner  and  outer  series  truncated  at 
300  are  nore  nearly  comparable. 

A  major  part  of  this  report  was  devoted  to  the 
development  of  ellipsoidal  harmonic  series,  in  particular, 
the  transformation  between  ellipsoidal  and  spherical 
harmonic  coefficients.  Although  the  downward  continuation 
error  in  ellipsoidal  series  is  generally  less  than  in 
spherical  series,  especially  in  the  polar  regions,  there 
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seeas  to  be  no  need  in  practice  (i<  300)  to  sake  the 
conversion  for  siaple  evaluations.  On  the  other  hand, 
because  the  spectral  coaponents  of  the  potential  on  the 
earth's  surface  bear  a  closer  reseablance  to  ellipsoidal 
haraonics  than  to  spherical  haraonics,  the  analysis  of 
terrestrial  data  (including  altiaetry)  is  sore  correctly 
coapared  to  (or  coabined  with)  the  ellipsoidal  spectrua. 

For  exaaple,  the  analysis  of  a  global  set  of  geoid 
undulations  (in  the  ellipsoidal  coordinate  systea  (&,X  )) 
yields  haraonic  coefficients  which  should  be  transf ora^d 
according  to  (3.34) ,  (3.35)  before  coaparing  thea  to 
potential  coefficients  derived  froa  satellite  data. 

The  expansions  of  the  sianlated  surface  and  gravity 
field  were  restricted  to  teras  of  degree  no  greater  than  300 
because  of  liaited  coaputer  storage  capabilities,  rather 
than  a  concern  about  excessive  coaputer  tiae.  Obviously, 
for  higher  expansions  of  the  potential,  the  error  analysis 
aust  be  redone  since  extrapolations  on  the  basis  of  Table  5 
are  risky.  In  any  new  study  of  the  downward  continuation 
error  one  should  endeavor  to  devise  a  density  distribution 
(such  as  point  aasses)  for  which  the  potential  function  can 
be  evaluated  to  any  accuracy,  thus  allowing  a  aore 
definitive  assessaeat  of  the  series  divergence.  Should  the 
error  ever  prove  to  be  relatively  significant,  it  is 
recoaaended  that  corrections  not  be  applied  to  spherical 
haraonic  coefficients  if  the  conversion  to  ellipsoidal 
haraonics  eliainates  the  significance  of  the  error. 

The  investigations  in  sections  2  and  3  have  relied  on 
approximate  foraulas  for  the  gravity  anoaaly  and  geoid 
undulation  (or  height  anoaaly)  in  order  to  sinplify  their 
functional  r el a tin ship  to  the  disturbing  potential.  In 
section  4,  corrections  to  these  approximations  were 
developed  with  the  preaise  that  the  disturbing  potential  is 
a  known  quantity  (e.g.  in  series  fom)  and  with  a  relative 
accuracy  on  the  order  of  the  square  of  the  earth's 
flattening.  These  corrections  should  be  applied  to  ng,  ^ 
whether  they  are  evaluated  using  the  spherical  series 
(2.50) ,  (2.58)  or  the  ellipsoidal  series  (3.51) ,  (3.53) 

(taking  due  account  of  the  coordinate  systeas  involved) . 
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/  >PEMDIXES 


>nverqenoe  and  Divergence 


The  "proof  of  convergence"  by  Arnold  (1978)  can  be 
outlined  as  follows.  It  is  well  known  that  any  reasonably 
well  behaved  function  (not  necessarily  continuous)  on  a 
surface  is  expandable  as  a  series  of  surface  spherical 
ha ra on ic  functions.  Instead,  however,  one  may  ask  whether 
the  set  of  solid  spherical  harmonics 
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considered  as  surface  functions  of  two  variables  0,X  (i.e. 

r=r  (  9 , X  ) )  is  also  a  complete  set  for  functions  defined  on 
the  surface.  In  this  case,  the  series  for  the  surface 
potential 
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would  be  a  unifornly  convergent  series  everywhere  on  the 
surface.  Consider  now  the  aonotonically  decreasing  sequence 
r"‘M/r,n* ' ,  r* >  r,  which  is,  noreover,  bounded.  Hence 
Abel*s  convergence  criterion  can  be  applied  to  claia  the 
convergence  of 
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for  every  r*>r.  Since  the  function  V  (r*,  0  ,  X  )  thus 
defined  is  haraonic  and  satisfies  the  boundary  values,  by 
the  uniqueness  of  the  boundary-value  problea,  V  must  be 
earth's  potential.  Purthsraore,  by  the  uniqueness  of  the 
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spherical  hanooic  expansion  of  the  potential,  the  above 
series  coincides  vith  the  series  that  would  be  deterained  on 
the  bounding  sphere.  Thus,  to  prove  convergence  of  the 
spherical  harmonic  series  everywhere  on  and  above  the 
earth's  surface,  indeed  everywhere  above  the  Bjerhammar 
sphere,  it  is  enough  to  prove  the  completeness  of  the 
functions  (A.1) .  Arnold  failed  here  as  his  "proofs"  of 
completeness  rely  on  the  assumed  troth  of  the  conjecture. 

In  a  more  recent  paper  (Arnold,  1980) ,  the  completeness  of 
the  functions  (A.1)  is  supposed  to  be  proved  by  showing  that 
there  exists  no  function  p  such  that 


//  p2  do  >  0  (A. 4) 

o 


and  such  that 


//  P  zxaa  do  "  0  ,  all  n,m  (A. 5) 


But  this  is  a  necessary  and  sufficient  condition  for 
completeness  of  orthogonal  functions;  the  ZHM>  are  clearly 
not  orthogonal  and  the  proof  fails  again. 

A  "proof  of  divergence"  of  the  potential  series  at  the 
earth's  surface  was  presented  by  Horrison  (1970)  under  the 
assumption  that  the  zonal  coefficients  do  not  decay,  in 
magnitude,  faster  than  saae  fixed  negative  power  of  the 
degree  n.  However,  the  "proof"  also  relies  on  the  erroneous 
statement  that  the  spherical  harmonic  series  diverges  if  its 
subseries  of  zonals  diverges  (the  latter  was  shown  for  the 
assumed  coefficient  decay) .  The  fallacy  of  this  argument  is 
easily  demonstrated  by  the  example  of  the  alternating  series 
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which  converges  (to  the  value  In 2)  ,  but  the  nubseries  of 
even  terns  only,  or  odd  terns  only,  by  itself  diverges. 
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£•  Series 


£2£  Sgmtgfl  Ellipsoid 


Por  points  P  below  the  bound in g  sphere,  the  serrated, 
homogeneous  ellipsoid  is  decoaposed  into  regions  for  which 
r?  >  rsi  and  r,  <  r,;  -  substituting  the  series  (2.7)  and 
(2.8)  into  (2.24),  the  potential  is  then  expressible  as  a 
convergent  series: 
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where  r$  =sin  (rr  ,r4j  )  ;  if  rf  >  t,;  ,  then  the  second  integral 
vith  respect,  to  r  vanishes#  Equation  (B.1)  readily  converts 
to 
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where  the  eguatorial  syaaetry  has  been  invoked  (only  even 
zona  Is  appear  in  the  series)  and  cti  <  r,  <  r1A  ;  if  r.  >  r,. 
for  all  if  then  1  •  Performing  the  integrations  with' 
respect  to  r  and  using 
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we  obtain 
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where  (using  equation  (2.59)  of  (Hobson,  1965,  p.33)) 
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=  <  *  pn-l(cos  0i+i>  *n^1 

Pi(cos  9^  -  Pi  (cos  8i+1>  ,  n  -  0 

As  reference  potential,  we  aay  use  the  zero  and  second 
degree  teres: 
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Then  the  disturbing  potential  becoaes  for  points  P  outside 
the  bounding  sphere: 
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and  for  points  P  below  the  bounding  sphere: 
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n^2 
n  even 


Pn(cos0p)  + 


(B.8) 


N-l  r  4  -  N-l  .  _  ,  r 

"  l  (“K*P0i  +  T  l  + 

i-j+1  P  °  5  i=j+l  5  rp  5  rP 

*  P2(cos9p)) 


□sing  the  definition  of  the  gravity  anoaaly  (2.42)  ,  the 
corresponding  series  ace  readily  found  to  be 


4gP  *  4"K*rP  al4  <»+SHin+l>tJ0  (^>n+3 


n  even 


(B.9) 


•pn<cos9p)  »  rp  >  R 


and  for  rp  <  B 


4*p  -  4"^rp  (X  i&rmzm  <J0 

n  £ven 


A  "  1  r  NrX  ,n+2  ,r0  ^n-2  4<2n+l)  *  ^  , 

n-0  1Sm  i-W5^  ^si  ~  <n+SXn--a>  >  ]  * 


n#2 
n  even 


•P„(cos9  )  + 
n  p 


(B.  10) 


^lipa<cosep)) 
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